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Chapter  1 


INTRODUCTION 


1.1  Background  and  Motivation 

Two  dimensional  complex  potential  methods  have  been  used  for  some  time 
for  the  prediction  of  the  stress  field  near  the  tips  of  a  crack  in  an  elastic  medium. 
These  methods  have  proven  to  be  very  accurate  and  generally  lead  to  the  solution  of 
two  dimensional  boundary  value  problems  with  Cauchy  type  singularities.  Using  these 
methods,  the  study  of  the  interaction  between  inclusions  and  cracks  has  been  active  for 
many  years.  Applications  of  the  results  of  theoretical  analyses  of  cracks  interacting 
with  inhomogeneities  have  had  wide  effect  in  the  implementation  of  a  broad  range  of 
new  and  emerging  material  systems.  The  understanding  gained  by  the  solution  of 
appropriate  elasticity  problems  provides  researchers  with  insight  into  the  mechanisms 
of  strengthening  and  toughening,  as  well  as  material  damage,  because  of  the  presence 
of  material  defects.  The  mechanisms  of  crack  growth  in  composite  materials  and 
ceramics,  and  that  of  strain  hardening  in  metal  alloys,  are  primary  examples  of  the 
direct  application  of  the  results  of  the  study  of  crack  inclusion  interaction. 
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crack  near  an  inclusion  provides  the  ability  to  quantitatively  assess  the  energy  absorbed 
by  crack  path  deflection.  This  provides  a  basis  for  a  fundamental  elasticity  model  for 
material  toughening.  Having  such  a  method  will  provide  researchers  with  a  useful  tool 
to  investigate  the  sources  of  toughening  in  materials,  or  to  guide  material  scientists  in 
their  efforts  to  improve  the  fracture  toughness  of  present  or  new  materials.  It  is  just 
such  a  method  which  is  described  in  this  work,  with  appropriate  examples,  test  cases, 
and  verification  of  results  m  order  to  convince  the  reader  that  the  proposed  technique 
is  a  viable  one,  and  a  useful  one. 

Methods  to  predict  the  path  of  a  crack  have,  in  the  past,  been  based  upon 
the  Finite  Element  Method  using  specialized  crack  tip  elements  which  are  intended  to 
capture  the  singularity  at  the  crack  tip.  These  predictions,  by  their  nature,  have  been 
computationally  intensive.  At  each  increment  of  crack  growth,  a  new  finite  element 
mesh  must  be  generated.  The  crack  has  propagated  into  the  mesh  during  this 
increment  of  growth,  and  the  crack  tip  element  must  stay  at  the  crack  tip.  Therefore, 
at  each  increment  of  crack  growth,  significant  computation  has  to  be  undertaken  just 
to  remesh  the  vicinity  of  the  crack  tip  to  model  that  increment  of  growth. 

The  proposed  method  for  crack  path  prediction  is  based  upon  a  Boundary 
Integral  approach,  similar  to  a  reasonably  elementary  Boundary  Element  solution, 
where  the  boundary  in  question  is  the  crack  itself.  This  method  has  been  used 
successfully  with  two  dimensional  potential  theory  to  calculate  the  stress  field  around 
the  tip  of  the  crack,  and  therefore  can  be  used  to  predict  the  direction  of  crack  growth. 
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Once  the  direction  of  crack  growth  is  known,  the  crack  can  be  grown  by  an  increment 
in  the  specified  direction  very  simply  by  allowing  the  boundary  (the  crack)  to  grow  by 
that  increment  in  that  direction.  The  re-meshing  of  the  boundary,  then,  is  simply  the 
extension  of  a  line  by  a  known  increment,  in  a  predicted  direction.  If  the  crack  is 
parameterized  as  a  smooth  set  of  cubic  interpolating  functions  commonly  referred  to 
as  a  cubic  spline,  the  extension  can  take  on  any  orientation  required  to  meet  a  specific 
crack  growth  criterion.  This  parameterization,  using  cubic  splines,  and  a  solution  to 
account  for  the  fact  that  the  naturally  growing  crack  is  not  straight,  are  developed  in 
detail  in  this  work. 

One  of  the  major  side  benefits  of  this  type  of  prediction  model  is, 
therefore,  that  it  is  much  less  computationally  intensive,  and  can  be  used  as  the  basis 
for  performing  sensitivity  studies,  or  for  a  stochastic  model  of  material  toughening. 
Of  course,  the  major  attraction  of  this  method  is  that  it  uses  an  exact  elasticity  solution 
which  can  be  used  as  a  basis  for  a  fundamental  mechanics  model  of  material 
toughening. 


1.2  Survey  of  Past  and  Current  Literature 

The  field  of  Linear  Elastic  Fracture  Mechanics  has  had  a  long  and 
venerable  tradition.  Early  work  by  Griffith  [1],  and  Irwin  [2]  laid  the  fundamentals 
for  what  has  become  an  indispensable  tool  for  the  practicing  Mechanical  Engineer. 
As  the  technique  has  become  established  and  accepted,  many  elasticians  have  solved 


3 


many  problems  of  interest.  Elasticity  solutions  have  been  mainly  for  two  dimensional 
regions,  both  for  their  simplicity,  and  for  the  ability  to  apply  two  dimensional  stress 
potentials  to  those  regions.  The  potential  methods  of  Muskhelishvili  [3]  have  been 
used  with  tremendous  success,  to  solve  some  very  complicated  and  difficult  elasticity 
problems  in  the  field.  Paris  [4]  and  Erdogan  [4,5],  among  others,  have  been 
instrumental  in  the  application  of  these  methods  to  solve  problems  of  significant 
interest  in  the  engineering  community.  These  problems  have  involved  interfaces, 
interactions,  welded  structures,  open  holes,  and  a  multitude  of  stress  problems  where 
the  stresses  become  singular  at  some  point  of  interest  in  the  field  in  question.  This 
work  deals  with  the  solution  of  one  of  these  problems  in  which  the  crack  is  interacting 
with  an  elliptical  inclusion,  and  is  growing  along  a  path  that  is  affected  by  the 
presence  of  that  inclusion. 

The  problem  of  the  interaction  between  a  circular  inclusion  and  a  static 
crack  has  been  solved  by  Atkinson  [6]  for  a  radial  crack,  and  by  Erdogan,  Gupta,  and 
Ratwani  [7]  for  an  arbitrarily  oriented  crack.  Erdogan  and  Gupta  [8]  later  solved  the 
problem  in  which  the  crack  crosses  the  interface.  These  solutions  are  based  on  the 
Green’s  Functions  derived  from  one  solution  of  a  dislocation  interacting  with  a  circular 
inclusion  (Dundurs  and  Mura  [9]  and  Dundurs  and  Sendeckyj  [10]).  Grief  and  Sanders 
[11]  used  the  same  techniques  to  solve  the  problem  of  the  effect  of  a  stringer  on  a 
cracked  sheet  of  material.  Their  solution  with  respect  to  this  work  represents  the  other 
extreme  in  ellipticity  ratio  from  that  done  by  Atkinson  for  the  circular  inclusion. 
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Santare  and  Keer  [12],  presented  the  two-dimensional  solution  for  the  dislocation 
outside  a  rigid  elliptical  inclusion,  using  the  complex  potential  methods  of 
Muskhelishvili  [3],  Their  Green’s  Function,  when  applied  to  a  straight  crack  using  the 
methods  described  in  this  work  matches  Atkinson’s  [6]  results  for  a  circular  inclusion, 
and  that  of  Grief  and  Sanders  [11]  for  the  line  inclusion  or  stringer.  Special  attention 
was  paid  in  [12]  to  the  rotation  of  the  elliptical  inclusion  and  the  effect  that  has  on  the 
stress  field  around  the  ellipse.  Results  were  compared  with  a  power  series  solution 
found  by  Stagni  and  Lizzio  [13],  which  did  not  take  into  account  the  rotation  of  the 
ellipse.  The  comparison  showed  that  the  rotation  has  a  significant  effect  on  the  stress 
field  in  many  instances.  In  another  paper,  Santare,  Keer,  and  Lewis  [14]  solved  a 
related  problem  of  an  elliptical  hole  at  a  bone/implant  interface,  with  symmetrical 
cracks  radiating  from  the  edges  of  the  ellipse,  along  the  x-axis.  Solutions  to  the 
resulting  singular  integral  equations  were  found  using  a  numerical  scheme  proposed 
by  Gerasoulis  [15].  The  technique  described  in  this  work  uses  the  potentials  calculated 
by  Santare  and  Keer  [12],  and  the  numerical  technique  described  by  Gerasoulis  [15] 
to  calculate  the  stress  intensity  factors  for  the  advancing  crack.  These  static  results 
were  reported  by  Patton  and  Santare  [16]  in  1990. 

Besides  the  interaction  problem,  the  stress  due  to  a  non-straight  crack  path 
is  also  important,  and  must  be  taken  into  account.  The  kinked  crack  problem  was 
attempted  by  several  researchers  in  the  early  1970’s,  with  some  success.  Palinaswamy 
and  Knauss  [17]  besides  providing  a  very  good  review  on  the  subject,  solved  the 
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problem  of  an  infinitesimal  branch  or  kink  off  the  tip  of  a  crack  using  a  Fourier  series 
approximation.  Their  results  are  in  agreement  with  an  earlier  solution  by  Bilby  and 
Cardew  [18],  and  have  been  found  by  others  to  be  correct.  Other  researchers 
(Dudukalenko  and  Romalis  [19],  and  Hussain,  Pu,  and  Underwood  [20])  came  up  with 
analytical  expressions  for  the  complex  potentials,  but  their  results  have  been  since 
found  to  be  inconsistent  with  the  results  published  by  other  authors.  Chateijee  [21], 
Gupta  [22],  and  Kitigawa,  Yukki,  and  Ohira  [23],  all  independently  solved  the  problem 
of  a  finite  extension  to  the  crack,  and  their  results  agree  with  one  another.  K.  K.  Lo 
[24]  finally,  in  a  paper  in  1978  presented  a  unified  approach  to  solving  the  kinked 
crack  problem  that  is  most  commonly  regarded  as  correct  and  most  commonly 
referenced  at  present.  His  results  match  those  of  [17]  and  [18]  for  the  infinitesimal 
kink,  and  the  results  of  [21],  [22]  and  [23]  for  the  finite  length  crack  extension. 

In  the  early  1970’s,  a  group  of  Russian  mathematicians  solved  the  problem 
of  a  crack  with  a  curved  extension  by  the  use  of  a  first  order  perturbation  method. 
Banichuk  [25]  found  the  first  order  perturbation  solution  and  came  up  with  complex 
potentials  suitable  to  a  solution  using  Muskhelishvili’s  [3]  method.  Goldstein  and 
Salganik  [26]  then  applied  this  set  of  potentials  to  the  solution  of  a  finite  crack  with 
a  curved  extension.  The  same  authors  reported  this  work  in  the  western  literature  [27] 
in  1974.  Since  that  time,  Cotterell  and  Rice  [28]  used  this  technique  to  solve  a  similar 
problem,  coming  up  with  a  somewhat  different  way  of  describing  the  problem.  They 
explicitly  calculate  the  stress  intensity  factors  without  the  need  to  write  integral 
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equations,  as  was  done  by  Goldstein  and  Salganik  in  [26],  Their  results  compare 
favorably  with  that  of  Goldstein  and  Salganik  [26].  Karihaloo  et  al.  [29]  included 
second  order  terms  in  the  same  solution  to  refine  the  prediction,  and  came  up  with  an 
elaborate  expression  for  the  stress  intensity  factors.  From  these  expressions,  the 
authors  derive  approximate  expressions  for  these  stress  intensity  factors  in  terms  of  an 
infinite  series.  Their  conclusions  from  this  work  were  that  the  path  of  the  propagating 
crack  depends  not  only  on  the  in-plane  normal  stress  as  concluded  by  Cottrell  and  Rice 
[28],  but  also  on  the  derivatives  of  the  stress  intensity  factors  with  respect  to  the 
original  crack  length.  They  go  on  to  conclude  that  a  crack  without  an  initial  kink  in 
a  non-homogeneous  stress  field  can  have  a  smooth  curved  path,  as  the  Mode  II  stress 
intensity  factor  will  be  changing  continuously  along  the  path  of  crack  growth,  under 
the  influence  of  the  non-homogeneous  stress  field.  That  is  exactly  the  situation  that 
is  examined  in  this  work,  that  of  a  crack  growing  through  the  inhomogeneous  stress 
field  caused  by  the  presence  of  a  rigid  elliptical  inclusion. 

Other  solutions  to  the  curved  and/or  kinked  crack  problem  have  been 
proposed  in  the  literature  in  recent  years.  Sumi,  Nemat-Nasser,  and  Keer  [30]  resolved 
the  first  order  perturbation  problem  originally  solved  by  Banichuk  [25]  and  Goldstein 
and  Salganik  [26],  with  the  additional  complexity  of  taking  into  account  the  boundaries 
of  the  region.  The  original  work  by  Banichuk  [25],  and  Goldstein  and  Salganik  [26] 
had  been  done  for  a  finite  crack  in  an  infinite  domain,  with  no  geometry  effects  of  the 
domain  taken  into  account.  Sumi,  Nemat-Nasser,  and  Keer  [31]  went  on  to  create  a 
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combined  analytical  and  finite  element  solution  to  the  crack  path  prediction  problem. 
Sumi  [32]  and  Sumi,  Ohashi,  and  Emura  [33]  have  gone  on  to  predict  crack  paths  in 
welded  structures  [32]  and  near  open  holes  in  cracked  sheets  [33].  The  second  of 
these  papers  included  an  experimental  investigation  of  cracks  interacting  with  open 
circular  holes.  Their  predictions  come  reasonably  close  to  their  experimental  results, 
but  more  importantly  for  this  work,  it  provides  experimental  results  with  which  to 
verify  the  predictions  made  with  the  present  technique.  Another  very  interesting 
solution  to  the  curved  crack  problem  has  been  proposed  by  Sur  and  Altiero  [34]  in 
which  they  do  not  formulate  the  problem  in  terms  of  the  derivative  of  displacement 
or  Burger’s  vector,  but  rather  in  terms  of  the  crack  displacement  itself.  They  claim 
that  their  technique  avoids  the  singularities  at  the  tips  and  kinks  of  the  crack,  and  that 
their  equations  are  not  hypersingular,  as  other  attempts  to  do  this  have  been  (see 
Ioakimidis  [35]).  They  compare  their  results  to  those  of  Kitigawa  et  al.  [23]  for  the 
kinked  crack,  with  good  agreement. 

The  criterion  which  is  used  to  predict  the  direction  in  which  a  crack  will 
kink  has  also  been  a  rather  controversial  problem.  There  are  basically  two  different 
schools  of  thought  regarding  the  correct  criterion  to  use  for  prediction  of  the  path  of 
an  advancing  crack,  stress-based  methods  and  energy-based  methods.  Erdogan  and  Sih 
[36]  proposed  a  criterion  which  has  come  to  be  called  the  "Maximum  Normal  Stress 
Criterion,"  in  which  they  state  that  a  naturally  growing  crack  in  an  isotropic  material 
will  grow  in  the  direction  in  which  the  stress  normal  to  the  crack  propagation  direction 


8 


is  greatest.  This  criterion  implies  that  the  stress  parallel  to  that  direction  would  be  a 
minimum,  or  that  the  shear  stress  in  that  direction  would  vanish,  as  the  crack  would 
propagate  in  a  direction  normal  to  the  maximum  principal  tensile  stress.  One  can 
easily  verify  this  by  rotating  a  plane  stress  field  into  its  principal  directions.  In  the 
principal  directions,  there  is  no  shear  stress.  The  authors  substantiate  their  claim  with 
experimental  data  which  seems  to  validate  their  analytical  results.  Other  authors  have 
applied  their  criterion  to  differing  crack  models.  McClintock  [37]  and  Cotterell  [38] 
noticed  a  discrepancy  when  the  elliptical  edge  crack  model  was  used  rather  than  the 
slit  model  used  by  Erdogan  and  Sih  [36].  Other  authors  (Williams  and  Ewing  [39], 
Finnie  and  Saith  [40],  and  Ewing  and  Williams  [41])  discussed  the  implication  of 
including  more  terms  in  the  near  tip  approximation  using  this  criterion.  Sih  [42] 
proposed  a  new  criterion  which  is  supposed  Lo  be  a  combination  of  a  stress  method 
and  an  energy  method,  which  he  called  the  Strain  Energy  Density  Criterion,  or  S- 
criterion.  This  criterion  has  come  under  some  criticism,  as  it  does  not  have  as  much 
fundamental  basis  as  the  energy  and  stress  methods. 

The  name.  Maximum  Energy  Release  Rate  Criterion,  was  initially  proposed 
by  C.  H.  Wu  [43],  in  1976  when  he  reviewed  all  of  these  works,  and  proposed  a 
criterion  for  fracture  which  is  consistent  with  Griffith’s  [1]  energy  release  concept. 
This  was  the  first  real  treatise  on  energy  based  methods  that  was  published.  This 
criterion  was  mentioned  by  Erdogan  and  Sih  [36],  and  a  solution  was  given  by 
Palinaswamy  and  Knauss  [17],  even  though  none  of  these  authors  explicitly  call  it  the 


Maximum  Energy  Release  Rate  Criterion.  In  his  work,  Wu  solved  the  problem  of 
describing  the  energy  released  by  an  infinitesimal  kink  at  the  end  of  a  straight  crack 
with  the  use  of  the  potential  methods  of  Muskhelishvili  [3],  and  then  applied  his 
criterion  that  the  energy  released  by  this  infinitesimal  propagation  should  be  a 
maximum  with  respect  to  the  direction  at  which  the  kink  propagates.  He  compares  his 
results  to  Erdogan  and  Sih  [36],  Sih  [42],  and  the  Fourier  series  solution  found  in 
Palaniswamy  and  Knauss  [17],  The  only  discrepancies  noted  are  for  a  crack  loaded 
in  pure  shear,  where  the  angle  of  the  kink  is  on  the  order  of  70  to  80  degrees.  In  the 
literature,  there  was  much  controversy  about  this  particular  problem,  but  in  the  work 
described  in  this  document,  the  kink  angles  are  very  small,  as  the  crack  is  growing  in 
a  stress  field  which  is  a  tensile  field  normal  to  the  original  crack.  The  fact  that  these 
two  major  criteria  produce  identical  results  for  small  kink  angles  in  a  stress  field  which 
does  not  produce  significant  Mode  II  loading,  leads  to  the  conclusion  that  either  is 
acceptable  for  the  work  described  in  the  following  pages,  and  the  major  reason  to 
choose  one  over  another  is  ease  of  implementation. 


1.3  Dissertation  Outline 

The  remainder  of  this  chapter  describes  the  content  and  organization  of  this 
dissertation,  in  enough  detail  that  the  reader  should  be  able  to  understand  the 
techniques  presented,  and  follow  the  methodology  presented. 

Chapter  2  presents  a  brief  summary  of  the  fundamentals  of  the  two 
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dimensional  potential  methods  used  in  deriving  the  crack  solution.  The  methods  of 
Muskhelishvili  [3]  are  presented  along  with  the  methods  by  which  the  singular  integral 
equations  used  to  calculate  stress  intensity  factors  are  derived  for  crack  solutions.  Also 
included  is  a  short  description  of  the  methods  of  solution  that  have  been  employed  in 
the  past  by  practitioners  of  Linear  Elastic  Fracture  Mechanics.  At  the  end  of  chapter 
2,  a  brief  discussion  of  toughening  mechanisms  in  materials  is  presented  in  order  to 
motivate  the  discussion  of  solutions  for  cracks  near  inhomogeneities. 

Chapter  3  presents  the  solution  to  the  point  dislocation  interacting  with  a 
rigid  elliptical  inclusion.  The  formulation  of  the  problem,  using  complex  potential 
methods,  is  presented,  and  the  potentials  that  were  calculated  by  Santare  and  Keer  [12] 
are  presented.  The  chapter  closes  with  a  summary  of  solution  technique  used  to  model 
a  straight  crack,  using  the  point  dislocation  potentials.  Results  of  that  analysis,  which 
are  original,  and  a  comparison  with  previous  work  are  presented.  This  chapter 
concludes  the  review  of  work  done  by  others,  with  the  discussion  of  the  solution  for 
a  dislocation  interacting  with  an  elliptical  inclusion.  The  solution  of  the  straight  crack 
problem,  and  all  subsequent  work  that  is  presented  in  chapters  4,  5,  and  6  is  original 
work. 

The  material  presented  in  chapter  4  is  intended  as  background  for  the 
technique  used  to  predict  the  path  of  the  crack.  The  most  fundamentally  difficult  pan 
of  this  is  the  ability  to  model  a  curvilinear  crack  in  a  complex  stress  field  (created  both 
by  the  applied  load  and  by  the  inclusion  itself).  The  chapter  begins  with  a  discussion 
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of  crack  kinking  and  curvature  solutions,  and  their  applicability  and  suitability  to  the 
chosen  method  of  solution.  Historically,  these  solutions  have  been  used  only  to 
correctly  calculate  the  stress  intensity  factors  for  a  non-straight  crack.  They  have  been 
intended  to  model  the  reaction  of  a  crack  to  a  change  in  the  direction  of  the  applied 
load.  There  has  been  little  done  to  extend  these  solutions  to  a  path  prediction  model, 
or  to  include  them  in  a  model  of  interaction  with  inhomogeneities.  The  chapter 
continues  with  the  description  of  a  perturbation  solution  which  was  performed  in  the 
early  1970’s  by  a  group  of  Russian  mathematicians  which  is  particularly  applicable  to 
the  problem  at  hand.  As  a  result  of  the  adoption  of  that  solution,  the  crack  must  be 
parameterized  in  such  a  fashion  as  to  yield  the  local  spatial  first  and  second  derivatives 
along  the  crack.  A  cubic  spline  parameterization  is  particularly  suitable  for  this 
parameterization,  and  that  is  also  described.  The  chapter  closes  with  a  description  of 
the  crack  extension  algorithm,  and  a  comparison  of  the  solution  with  known  kinked 
crack  solutions. 

Chapter  5  describes  the  method  used  to  predict  the  crack  path  near  the 
inclusion.  It  begins  with  a  historical  review  of  crack  propagation  laws,  and  a  rationale 
for  the  choice  of  the  Maximum  Normal  Stress  as  the  criterion  used  to  predict  the  local 
direction  of  incremental  crack  growth.  The  algorithm  used  to  predict  the  local 
direction  is  described,  and  a  verification  of  the  prediction  is  described.  The  chapter 
continues  with  parametric  studies  of  path  deflection  as  a  function  of  ellipticity  ratio  of 
the  ellipse,  and  on  the  initial  location  and  orientation  of  the  crack.  It  concludes  with 


12 


an  estimate  of  energy  absorbed  by  crack  path  deflection  for  several  crack  paths  of 
interest 

The  last  chapter  of  this  dissertation  reviews  significant  results  and  presents 
conclusions  obtained  from  the  study  of  crack  path  deflection,  along  with  a  few  ideas 
for  future  work  that  should  come  as  a  result  of  this  study. 
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Chapter  2 
FUNDAMENTALS 


2.1  Two  Dimensional  Complex  Potential  Methods 

It  can  be  shown  that  the  stress-strain  relations  for  two  dimensional  plane 
problems  can  be  written  in  indicial  notation  as  follows  (Muskhelishvili  [3]): 

ea6-— (o  B-x5  )  (2.1) 

where  Muskheli  lvili’s  constant,  K  =  (3-v)/(l+v)  lor  plane  stress,  and  K  =  3-4v  for 
plane  strain,  v  is  the  Poisson’s  ratio  of  the  material,  the  strain  is  defined  by  e^,  the 
stress  by  a^,  p  is  the  shear  modulus  of  the  material,  and  the  kronecker  delta  (5)  has 
its  usual  indicial  definition.  Following  the  complex  potential  methods  of 
Muskhelishvili,  it  is  assumed  that  the  complex  variable  z  defines  x  and  y  coordinates 
as  in  z  =  x  +  iy,  where  i  is  the  imaginary  number.  Furthermore,  displacements  and 
forces  can  take  on  complex  values  as  D  =  u  +  iv  where  u  is  displacement  is  the  x- 
direction  v  is  displacement  in  the  y-direction,  and  the  forces  as  F  =  F,  +  iFy,  where 
the  subscripts  refer  to  forces  in  those  directions.  From  these  quantities,  the  stresses 
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and  displacements  in  the  region  can  be  written  as  a  function  of  two  holomorphic 
potential  functions  4>  and  \j/  as  follows: 

°“+°»=2[|<l>(z)+|4,(z)1  (22) 

+  =2[Z“~~4>Cz)+  “■  tjf(z)]  (2.3) 

dz 

2p(w+iv)=K<i>(z)-z-7-4>(z)-'Kz)  (2-4) 

dz 

In  the  stress  analysis  of  a  two  dimensional  region,  therefore,  what  is 
required  is  the  calculation  of  these  two  potential  functions.  In  a  great  number  of  cases, 
these  two  potentials  can  be  calculated  by  evaluating  two  Cauchy  integrals,  as  given  by 
Muskhelishvili.  In  the  case  of  a  single  dislocation  in  an  infinite  homogeneous  elastic 
medium,  these  integrals  give  the  following  potentials: 

<|)d(z)=Ylog(z  -Zq)  (2  5) 

'Mz)=Yk>g(z-Z0)-Y-^-  (2'6) 

z_zo 

where  the  dislocation  is  included  in  the  form  y  =  p(bx+iby)/in(K+l),  z  is  some  point 
in  the  complex  plane  and  z,  is  the  position  of  the  dislocation.  The  dislocation,  in  the 
form  described  above  (b,  and  by),  is  called  the  Burger’s  vector,  which  is  related  to  the 
derivative  of  the  displacement  at  the  point  of  the  dislocation.  Obviously,  the  stresses 
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become  singular  as  one  approaches  the  dislocation,  as  would  be  expected.  Therefore, 
to  study  the  stress  field  near  a  dislocation,  or  along  the  line  of  dislocations  commonly 
used  to  model  a  crack,  one  must  deal  with  the  singularity  in  stress.  In  linear  elastic 
fracture  mechanics,  the  stresses  along  the  crack  are  assessed  by  modelling  the 
displacement  of  the  two  faces  of  the  crack  as  a  line  of  dislocations  or  crack  opening 
displacements  which  is  unknown.  The  stresses  which  are  related  to  the  strains,  or 
derivatives  of  displacement,  through  the  stress-strain  relations  (equation  2.1),  therefore, 
can  be  written  as  a  set  of  singular  integral  equations  where  the  unknown  distribution 
of  dislocations  is  inside  the  integrals.  The  integrals  are  taken  along  the  crack,  and  the 
value  of  the  integrands  are  singular  at  the  end  points,  namely  the  tips  of  the  crack. 


2.2  Solution  Methods 

Several  techniques  have  been  used  in  the  past  to  solve  the  singular  integral 
equations  that  result  from  the  stress  analysis  of  a  crack,  but  the  most  commonly  used 
is  an  approximation  of  the  singular  part  of  the  integral  with  a  suitable  polynomial 
(Chebychev  or  LaGrange  are  common  choices),  and  a  direct  numerical  evaluation  of 
the  non-singular  part  of  the  integral.  The  technique  developed  by  Gerasoulis  [15] 
which  is  noted  in  the  discussion  of  pertinent  literature  in  the  introduction  is  the  one 
employed  in  this  work. 

Since  the  stress  analysis  is  linear,  the  principle  of  superposition  applies. 
Therefore,  if  the  problem  in  question  can  be  modeled  as  a  superposition  of  several 
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stresses  upon  one  another,  the  stresses  due  to  each  of  them  can  be  calculated,  and  then 
added  up.  In  the  analysis  of  cracks  interacting  with  inclusions,  much  use  of  the 
principle  of  superposition  is  made.  In  practical  use,  the  methods  to  evaluate  the 
singular  integral  equations  have  a  singular  part  which  represents  the  crack  tip 
singularity,  and  a  non-singular  part  which  represents  the  superposition  of  all  of  the 
other  stresses  which  act  on  the  crack  at  some  distance  from  the  crack  tips.  Therefore, 
to  analyze  the  effect  of  an  inclusion  interacting  with  a  crack,  one  needs  to  calculate 
potentials  <J>  and  \j/  for  the  interaction  between  a  dislocation  and  an  inclusion.  Using 
these  potentials,  the  stresses  are  calculated,  superposed  and  summed  along  the  crack 
length,  resulting  in  an  integral  equation  valid  along  the  length  of  the  crack.  Therefore, 
the  integrals  typically  take  the  following  form  for  a  finite  crack: 

Cx  f-^-d^+Cz  jKizjJHzotdzo^Fiz)  (2-7) 

zi  z  zi 

where  the  two  constants  (C,  and  Q)  are  based  upon  material  and  geometric  properties, 
the  integrals  are  taken  between  the  two  crack  tips  z,  and  z^  and  the  dislocation 
density,  b(zo),  is  an  unknown  function  along  the  crack  which  is  evaluated  at  the  points 
Zq.  This  form  explicitly  separates  the  singular  part  of  the  integral  equahon  as  the  first 
integral,  and  the  non-singular  pan,  or  kernel  (KU.Zq)),  as  the  second  integral.  The 
forcing  function,  or  the  load  applied  to  the  crack,  is  given  as  F(z)  on  the  right  hand 
side  of  the  integral  equation. 
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2.3  Toughening  Mechanisms  in  Materials 

One  of  the  most  common  methods  of  strengthening  and  toughening 
engineering  materials  has  been  the  incorporation  of  inhomogeneities  into  the  base 
structure  of  the  material.  Composite  materials  are  probably  the  most  well  known 
example  of  the  purposeful  use  of  this  method  for  strengthening.  In  an  organic 
composite  material,  the  relatively  weak  and  brittle  organic  matrix  material  is  both 
strengthened  and  toughened  by  the  incorporation  of  fiber  reinforcement.  In  a  typical 
discontinuously  reinforced  metal  matrix  composite  material,  the  hard,  tough  reinforcing 
phase  is  interspersed  within  the  metallic  grain  structure.  The  same  is  true  of  ceramic 
matrix  composites,  where  a  harder  ceramic  reinforcement  is  dispersed  into  the  softer 
ceramic  matrix  surrounding  it.  A  more  mundane  example  of  this  toughening  is  that 
of  high  strength  steels  and  aircraft  grade  aluminum  alloys.  In  both  of  these  materials, 
hard  intermetallic  compounds  are  precipitated  out  into  the  boundaries  around  the  grains 
of  the  metal,  and  strength  is  increased.  In  many  cases,  if  the  bonding  between  the 
matrix  material  and  the  reinforcement  is  good,  and  if  the  reinforcement  takes  on 
certain  geometric  and  physical  properties,  the  toughness  of  the  material  is  increased. 

One  of  the  more  compelling  theories  that  has  been  proposed  to  explain  this 
phenomenon  is  that  the  deflection  of  the  path  of  a  crack  through  the  material  absorbs 
a  significant  amount  of  energy,  thereby  toughening  the  material.  Rubinstein  [45] 
discusses  this  mechanism  at  length,  giving  solutions  to  the  problem  of  calculating 
toughness  increase  due  to  crack  path  and  crack  shielding  effects.  This  work  is  used 
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in  a  later  section  of  this  document  to  estimate  the  energy  absorbed  by  cracks  which 
are  effected  by  elliptical  inclusions.  It  would  seem  that,  from  this  work  and  that  of 
Rubinstein,  a  fundamental  model  of  this  crack  path  deflection  that  would  allow  the 
estimation  of  energy  absorbed  by  crack  path  deflection  would  be  useful  in 
investigations  into  material  toughening.  It  is  precisely  this  mechanism  that  provided 
the  motivation  for  the  work  that  is  described  in  this  document. 
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Chapter  3 


ELLIPTICAL  INHOMOGENEITY 


3.1  Problem  Formulation 

The  problem  of  calculating  the  interaction  between  a  dislocation  and  an 
elliptical  inhomogeneity  can  be  solved  by  determining  suitable  potential  functions  <J) 
and  y  such  that  the  classical  methods  described  in  this  work  can  be  used.  These 
potential  functions  can  be  calculated  by  solving  an  appropriate  boundary  value  problem 
that  can  be  cast  as  two  integrals  in  the  complex  plane.  This  problem  was  solved  by 
Santare  and  Keer  [12],  and  the  salient  points  are  described  in  the  following. 

The  geometry  of  the  problem  is  simplified  by  uang  a  conformal  mapping 
of  the  region  outside  an  ellipse  onto  the  region  outside  the  unit  circle,  as  shown  in 
Figure  (3.1).  The  function  that  does  this  mapping  can  be  written  as, 

Z=co(0=fl(C+-y)  (3-1) 

where 
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(3.2) 


and  the  parameter  which  can  be  thought  of  as  the  ellipticity  of  the  ellipse  is  given  by 


m=- 


a-b 

a+b 


(3.3) 


where  a  and  b  are  the  semi-axes  of  the  ellipse  shown  in  Figure  (3.2).  This  figure 


Figure  3.1.  Mapping  function  for  the  dislocation-inclusion  interaction  problem. 


displays  the  geometry  of  the  problem  of  a  crack  interacting  with  an  elliptical  inclusion. 
In  Figure  (3.2),  the  two  crack  tips  are  shown  as  Z,  and  Zj,  the  distance  from  the  ellipse 
to  the  close  crack  tip  is  defined  as  d,  the  angle  that  the  crack  makes  with  respect  to 
the  x-axis  is  0,  and  the  angle  that  the  close  tip  of  the  crack  makes  with  respect  to  the 
major  axis  of  the  inclusion  is  a.  The  parameter  that  defines  ellipticity,  m,  takes  on  a 
value  of  zero  for  a  circle,  and  a  value  of  one  for  a  line  along  the  x-axis.  A  value  of 
m  of  -1  describes  a  line  along  the  y-axis. 
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The  solution  for  a  dislocation  interacting  with  a  rigid  elliptical  inclusion 
is  used  to  formulate  the  Green’s  functions  for  the  crack  problem.  The  potentials  were 
calculated  by  Santare  and  Keer  [12],  and  are  restated  here  for  completeness.  They 
include  the  dislocation  in  the  form,  described  previously,  of  y  =  p  (bx  +  ibv)/i7t  (k  + 
1),  the  mapped  coordinates  C,  and  which  are  the  transformed  z  and  z0,  and  the 
rotation  of  the  ellipse  eo. 


4»(C)  =  ^logWC-Co+m/t-zn/Co)] 

K 


,  r(C‘1/Co>i  ,  r 

+log[ - - — ]  -  ylog[- 

j_  (1/Cj-w/Cp)  (1/Cq-Cq) 

(VTo-Tjm)  (\/T0-O 


(C-m/Q 

c  J 


+ 


kC 


(3.4) 
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where 


iKO  =  Ylog[J?(C-Co+w/C-/n/Co)]  -  Y 


Co  +  ml  C0 
C-C0+m/C-m/C0 


r  (C-l/Co) ,  .(C-m/Co), 

KYlogI — - — ]  -  ylog[ - - — ] 


(m/C0-l/Co)  (m/C0-C(//n)  „  . 

+  -  +  2uien— 

(m/Co-Co)  (m/C0-C)  °C 

»  1  +m(2  ,  y_  1/Cp  _  m/C0 

c2-m  k  c(c-i/Co)  Y  CiUC-m/Q 

0/Cjf  >  (1/Co-Co) 

+  JL  _zzrJo - __  . 

™  (1/Co-Co/m)  (C-l/Co)2  K^2 


(3.5) 


B  .r  ky  —  m  m2 
e0  =  Re  i[-^=-+y— +Y— 


m 

-Y—  -  ym 

kC0 


C0  Co  '•o 
(m/C0-l/O  (m/Co-C^m) 


(m/Co-Co) 

7  (l/C^-m/Co)  (1/Co-Co) 


(1/Co- Co/m) 


]/2pl?(l+— ) 

K 


(3.6) 


The  first  term  in  0  (equation  3.4)  and  the  first  two  terms  in  y  (equation 
3.5)  represent  the  potentials  for  dislocation  in  an  unbounded  medium,  and  are  the 
restatement  of  equations  2  .5  and  2.6  in  the  mapped  plane.  These  terms,  when  applied 
to  equations  2.2  and  2.3  to  calculate  the  stresses,  become  singular,  and  constitute  the 
first  integral  in  equation  2.7.  The  remainder  of  the  terms  in  both  potentials  account 
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for  the  interaction  between  the  ellipse  and  the  dislocation.  These  terms,  when  used  in 
equations  2.2  and  2.3,  constitute  the  nonsingular  stresses  due  to  the  interaction  between 
the  dislocation  and  the  inclusion.  Note  also  that  the  rotation  of  the  ellipse  is  taken  into 
account  through  these  nonsingular  terms  and  is  called  E0.  The  last  term  in  the 
expression  for  \y  (the  bracketed  terms)  contains  the  first  derivative  of  the  interaction 
portion  of  the  first  potential  0.  This  is  a  natural  consequence  of  the  method  of 
Muskhelishvili  [3]. 
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3.2  Solution  Method  for  a  Straight  Crack 


Summing  these  stresses  along  the  crack  length  (Figure  (3.2))  and  setting 
them  equal  to  the  stresses  due  to  the  external  load,  two  singular  integral  equations 
result.  The  solution  of  these  equations  will  give  the  unknown  distribution  of 
dislocations  along  the  crack,  from  which  the  stress  intensity  factors  at  the  crack  tips 
can  be  directly  calculated.  If  the  stresses  are  resolved  into  components  normal  to  the 
crack  (Mode  I),  and  along  the  crack  (Mode  II),  the  two  equations  can  be  stated  as 
follows: 


*2  ,  ,  v  *2 

o  +  /  o)^o  =  W 

z,  z  *o 


hb(zJ  * 

f-^dz o  +  /  a  t(z^)bt(z<)dz^  =  ft(zj 

t  z_zo  t 


where  the  n  and  t  subscripts  refer  to  the  normal  and  tangential  components  with 
respect  to  the  crack,  f„  and  f,  are  the  stresses  due  to  the  external  load,  and  the  non¬ 
singular  parts  of  the  potentials  are  denoted  as  a„  and  a,.  The  first  term  in  each  of 
these  equations  represents  the  Cauchy  singular  portion  of  the  stresses  (as  noted  above), 
and  the  second  term  contains  the  nonsingular  parts.  The  normal  and  tangential 
cartesian  stress  components  can  be  determined  using  a  standard  trigonometric 
transformation.  The  solution  to  this  problem  is  non-unique,  however,  until  one  more 


condition  is  met,  that  of  crack  closure.  Simply  stated,  the  crack  must  close,  or  the 
crack  opening  displacements  must  be  identically  zero  at  the  two  endpoints  or  crack 
tips.  This  condition  can  be  stated  in  the  formulation  of  this  problem  as  follows: 

/tGoVfeo  -  0  0-9) 

*1 

If  this  condition  is  met,  there  will  be  a  single,  unique  solution  to  the  unknown 
distribution  of  dislocations. 

Furthermore,  if  equations  3.7  and  3.8  are  rewritten  in  terms  of  the  polar 
form  of  a  vector  from  z  to  Zq  as  follows: 

£  =  =  pe 16  (3.10) 

where  the  singularity  is  explicitly  shown  as  p  goes  to  zero,  the  stress  due  to  the 
Cauchy  singular  portion  of  the  above  integral  equation  can  be  determined,  in  terms  of 
the  angle  0  (the  angle  of  the  crack,  as  shown  in  figure  3.2),  as  follows: 

o  =  —  [b  (cos0  +  cos30)  - 
*  P  y  (3.11) 

bx( 3sin0  +  sin30)] 

o_  =  —  [fc(3cos0  -  cos0)  - 

9  y  (3.12) 

£x(sin30  -  sin0)] 
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(3.13) 


=  —  [£>y(sin30  -  sin0)  + 
bx(  cos30  +  ccs0)] 

This  permits  calculation  of  the  x  and  y  components  of  the  dislocation  vector  (Burger’s 
vector)  for  a  straight  crack  in  any  orientation. 

Once  these  integral  equations  have  been  written,  a  numerical  treatment 
must  be  used  to  solve  for  the  unknown  distributions  of  dislocations  in  the  x  and  y 
directions,  b,  and  by.  There  are  several  treatments  to  this  prob^m  that  have  been 
published  in  the  literature,  all  of  which  use  a  polynomial  appioximation  of  the 
unknown  function.  The  numerical  technique  used  in  this  work  was  proposed  by 
Gerasoulis  [15]  in  1980.  This  technique  uses  a  piecewise  quadratic  polynomial 
representation  of  the  singular  and  non-singular  parts  of  the  integral  equation,  and 
reduces  the  integrals  to  a  discretized  set  of  algebraic  equations  well  suited  to  a  matrix 
type  solution.  However,  in  order  to  use  the  technique,  first  some  assumptions  about 
the  nature  of  the  singularities  at  the  tips  of  the  crack  must  be  made.  Following  the 
method  of  fracture  mechanics,  it  can  be  shown  that  the  stresses  at  the  crack  tips 
become  singular  as  \/pm.  The  dislocation  density  will  therefore  also  exhibit  this 
square  root  singularity  at  both  ends,  and  can  be  written  in  the  following  form: 

b(s)  =  S®-  (3.14) 

\Jl-S2 
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where  s  =  2(z  -  z,)/^  -  z,)  -  1.  The  functions  g(s)  are  therefore  continuous  and 
bounded  along  the  interval  s  e  [-1,1].  Also,  in  general,  the  solution  to  the  singular 
part  of  the  integral  equation  is  not  unique,  and  an  additional  condition  is  required.  In 
the  present  work,  that  condition  is  that  the  crack  must  close  at  the  crack  tips. 
Mathematically,  this  means  that  the  value  of  the  integral  along  the  interval  from  -1  to 
1  of  the  unknown  distribution  of  dislocations  (g(s))  must  be  zero. 

For  a  given  crack  discretization  parameter  n,  the  numerical  technique 
provides  2n+l  integration  points  and  2n  collocation  points,  as  follows: 

2JV-1 

E  wy  -  *  Atp  (3 15) 

7-lA  ...Off 

where  the  w’s  are  the  weighting  functions  that  are  applied  to  the  singular  pan  of  the 
integral,  the  v’s  are  the  coefficients  of  the  non-singular  kernel  (K),  the  i’s  are  the 
integration  points,  and  the  j’s  are  the  collocation  points.  The  weight  functions  w,(s) 
and  Vj  are  calculated  by  using  the  LaGrange  inteipolation  formula  for  three  points  to 
approximate  the  unknown  bounded  function  in  the  case  of  the  w’s,  or  the  unknown 
function  multiplied  by  the  kernel  (K(s,t))  in  the  case  of  the  v’s,  with  the  singularity 
explicitly  removed,  and  integrating  (summing)  along  the  interval.  The  final 
expressions  for  the  coefficients  that  make  up  the  weight  functions  are  rather 
complicated,  and  are  found  in  the  paper  by  Gerasoulis,  and  therefore  will  not  be 
restated  here. 
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The  calculation  of  the  unknown  g(s),  then,  amounts  to  the  creation  of  a 
matrix  of  the  weighting  functions,  w(s),  added  to  the  non-singular  kernel  multiplied 
by  the  weighting  functions  v,  in  the  form  [A] { x }  =  { b} .  As  stated  previously,  the 
matrix  [A]  is  not  a  square  matrix,  or  there  are  2n+l  equations  with  2n+2  unknowns 
until  the  condition  of  crack  closure  is  added  by  multiplying  the  non-singular  weight 
function  by  the  cosine  or  sine  of  the  local  angle  and  placing  that  value  in  the  last  row 
of  the  matrix,  as  follows: 

A(2N+2  ,y)=ki(/)  C„  (3.16) 

where  Cn  denotes  the  sine  or  cosine  of  the  local  angle  of  the  crack.  The  solution  then 
requires  inversion  of  the  matrix  [A],  and  multiplying  by  {b},  the  forcing  function,  to 
determine  the  unknown  { x }  or  in  this  case,  { g(t) } .  Gerasoulis  states  that  the  method 
works  equally  well  for  unequal  meshes,  but  in  this  work,  the  spacing  was  kept  equal 
for  simplicity  in  coding  of  the  algorithm.  A  value  of  n  above  6  was  found  to  be 
unnecessary  for  crack  lengths  on  the  order  of  the  size  of  the  inclusion. 


3.3  Results  and  Comparison  to  Previous  Work 

Figures  (3.3)  through  (3.9)  are  shown  to  compare  and  contrast  the  results 
of  the  solution  of  a  straight  crack  interacting  with  an  elliptical  inhomogeneity  with 
those  previously  published.  These  results  are  original,  and  are  reported  in  Patton  and 
Santare  [16].  In  all  of  these  figures,  the  results  are  presented  in  terms  of  the 
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normalized  Mode  I  stress  intensity  factor.  This  stress  intensity  factor  is  normalized  as 
follows: 

K, 

normalized  k.  =  -  (3.17) 

o  JL 

where  the  stress,  is  the  applied  load  on  the  infinite  medium,  the  crack  length  is  L, 
and  the  Mode  I  stress  intensity  factor  is  in  the  numerator.  The  normalization  in  this 
equation  means  that  without  the  presence  of  an  inclusion,  the  normalized  stress 
intensity  factor  will  have  a  value  of  one,  or  the  stress  intensity  factor  will  be  the  same 
as  that  for  a  crack  in  an  infinite  medium  without  an  inclusion. 

Since  the  greatest  interest  in  this  work  is  with  the  crack  tip  closest  to  the 
inclusion,  all  results  reported  are  for  the  stress  intensity  factor  for  the  tip  closest  to  the 
inclusion.  Also,  to  make  comparisons  between  different  results  meaningful,  some  of 
the  specific  geometry  parameters  are  held  fixed  for  the  reported  results.  Specifically, 
the  value  of  R  (the  size  of  the  inclusion)  was  held  fixed  at  unity,  as  was  the  length  of 
the  crack.  For  results  reported  for  a  rigid  inclusion,  Muskhelishvili’s  constant,  K,  is 
held  at  a  value  of  1.67.  Also,  the  shear  modulus,  p,  does  not  require  formal 
specification  as  a  result  of  the  normalization  procedure.  Therefore,  as  m  (the  ellipticity 
ratio)  is  varied,  the  values  of  a  and  b  from  Figure  3.2  can  be  calculated  using 
equations  3.2  and  3.3  for  all  of  the  results  presented  in  this  work.  It  was  discovered 
during  the  course  of  this  investigation  that  a  crack  length  of  approximately  the  size  of 
the  inclusion  provided  results  both  to  compare  with  previous  results,  and  also  to  make 
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the  effect  of  the  other  crack  tip  on  the  results  nearly  negligible.  By  this  it  is  meant 
that  the  results  for  a  unit  length  crack  interacting  with  a  unit  circular  inclusion  gave 
results  that  varied  very  little  from  a  crack  of  length  10  or  100,  when  normalized 
according  to  equation  3.17,  interacting  with  the  same  inclusion.  Since  the  intention 
is  finally  to  create  a  path  prediction  model,  with  the  crack  tip  clos  ;st  to  the  inclusion 
exhibiting  growth,  this  result  is  meaningful  in  that  the  effect  of  the  far  crack  tip  is 
minimized  in  that  prediction  model. 

The  problem  of  a  circular  inclusion  was  treated  by  Atkinson  [6]  and  by 
Erdogan  and  Gupta  [8],  and  that  for  the  elastic  line  inclusion  by  Grief  and  Sanders 
[11].  These  are,  of  course,  both  extremes  of  the  solution  for  the  elliptical  inclusion, 
and  provide  a  good  comparison.  In  the  paper  by  Grief  and  Sanders,  the  integral 
equation  is  formulated  somewhat  differently,  in  that  only  the  interaction  terms  are 
taken  into  account.  The  present  solution  takes  into  account  not  only  the  interaction 
terms,  but  also  the  stress  singularities  at  the  tips  of  the  crack.  The  results  of  that  paper 
therefore  have  to  be  added  onto  the  stress  due  to  the  singularities  to  provide  a 
comparison.  That  has  been  done,  and  the  results  for  a  rigid  inclusion  with  a  range  of 
ellipticities  (m)  are  shown  in  Figure  (3.3).  In  this  figure,  the  crack  is  oriented 
horizontally,  aligned  along  the  real  axis  (x-axis),  and  is  at  varying  distances  (d)  from 
the  edge  of  the  inclusion.  The  curve  in  the  figure  labeled  as  m=0  (a  circular  inclusion) 
reproduces  the  results  of  Atkinson.  A  comparison  of  the  two  sets  of  results  is  given 
in  Table  3.1.  In  the  work  by  Grief  and  Sanders,  rather  than  presenting  stress  intensity 
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Figure  3.3.  Stress  intensity  factor  versus  distance  from  ellipse  for  a  radially  oriented 
crack  and  varying  ellipticity  ratios  (load  normal  to  crack,  R=l.,  L/R=l„  k=1.67). 


factors,  they  present  a  ratio  of  the  stress  if  the  interaction  between  the  crack  and  the 
line  inclusion  (stringer)  is  taken  into  account  versus  the  stress  if  that  interaction  is  not 
present.  They  also  present  results  for  an  elastic  stringer,  rather  than  a  rigid  one. 
Therefore,  if  in  the  present  analysis,  such  a  ratio  is  calculated,  it  would  represent  the 
limit  as  their  analysis  approached  a  rigid  stringer.  Those  results  are  shown  in  Table 
3.2,  where  their  flexibility  parameter,  A  is  varied  from  a  value  of  10  (flexible  stringer) 
to  0. 1  (stiff  stringer),  and  compared  to  the  rigid  or  limiting  case  presented  in  this  work 
in  which  it  would  have  the  value  zero.  These  results  are  reasonable  and  follow  the 
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Table  3.1.  Comparison  of  present  results  with  those  of  Atkinson 


Distance  From  Inclusion 

Atkinson 

Present  Solution 

.1 

.43 

.4276 

.2 

.58 

.5867 

.3 

.68 

.6857 

.4 

.76 

.7527 

.5 

.80 

.8002 

.6 

.84 

.8353 

.7 

.86 

.8619 

.8 

.88 

.8825 

.9 

.90 

.8988 

1.0 

.91 

.9119 

trend  reported  in  the  Grief  and  Sanders  paper. 

Figure  (3.3)  clearly  shows  the  effect  on  the  stress  intensity  factor  due  to 
the  presence  of  the  inclusion,  and  the  effect  of  both  the  shape  and  orientation  of  the 
inclusion.  For  those  inclusions  that  are  long  and  slender,  aligned  normal  to  the  crack 
(negative  m’s),  the  effect  called  crack  shielding  is  predominant.  This  phenomenon,  the 
lessening  of  the  Mode  I  stress  intensity  factor  due  to  the  presence  of  an  inclusion, 
called  crack  shielding,  is  exactly  the  phenomenon  postulated  to  be  the  toughening 
mechanism  attributed  to  inclusions.  It  is  also  interesting  to  note  that  there  is  slight 
anti-shielding,  or  a  slight  increase  on  K,  when  a  line  inclusion  is  oriented  in  the  same 
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Table  3.2.  Comparison  of  present  results  with  those  of  Grief  and  Sanders 


Distance  From 
Inclusion 

Grief  and 

Sanders  (A  = 

10.) 

Grief  and 
Sanders  (A  = 

•  1) 

Present  Work 
(A  =  0) 

.05 

.92 

.67 

.5279 

.1 

.96 

.78 

.6619 

.15 

.975 

.85 

.7517 

.2 

.98 

.88 

.8159 

.25 

.985 

.92 

.8627 

.3 

.99 

.94 

.8974 

.35 

.99 

.955 

.9232 

4. 

.99 

.96 

.9428 

direction  as  the  crack.  In  a  path  prediction,  one  would  expect  that  the  rounder 
inclusions  or  inclusions  oriented  normal  to  the  path  of  the  crack  would  deflect  the 
crack,  whereas  a  line  inclusion,  or  a  long,  thin  inclusion  oriented  in  the  same  direction 
as  the  crack  would  tend  to  attract  the  crack,  because  of  this  shielding  or  anti-shielding 
effect. 

Figure  (3.4)  shows  similar  results  to  Figure  (3.3)  except  that  in  this  Figure, 
the  inclusion  is  an  open  hole.  Muskhelishvili  [3]  noted,  in  his  work,  that  if  the 
material  parameter,  K,  was  set  to  a  value  of  -1,  which  is  physically  unrealistic,  the 
solution  to  the  rigid  inclusion  problem  becomes  that  for  an  open  hole.  In  Figure  (3.4), 
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Figure  3.4.  Stress  intensity  factor  versus  distance  from  open  hole  for  a  radial  crack 
and  varying  ellipticity  ratios  (load  normal  to  crack,  R=l.,  L/R=l.,  k=-1.) 

then,  the  parameter,  K,  is  set  to  a  value  of  -1,  and  the  results  are  for  open  holes  with 
varying  ellipticity  ratios.  Note  that  the  character  of  this  figure  is  opposite  that  of  the 
previous  figure,  in  which  the  inclusion  is  rigid,  in  that  as  the  close  crack  tip 
approaches  the  open  hole,  the  stress  intensity  factor  increases.  If,  therefore,  an 
inclusion  would  be  expected  to  deflect  a  crack  because  of  shielding,  a  hole  would  be 
expected  to  attract  a  crack  because  of  anti-shielding.  This  is  physically  reasonable,  and 
has  been  seen  in  practice. 

Figures  (3.5)  through  (3.9)  show  interesting  results  of  the  solution  for 
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differing  cracks  at  differing  orientations  interacting  with  differing  inclusions.  Figure 
(3.5)  shows  results  for  a  radial  crack  (0  =  a  in  Figure  3.2)  at  several  different 
orientations  interacting  with  an  inclusion  with  varying  ellipticity  ratios.  The  distance 
from  crack  to  inclusion  is  held  constant  in  this  Figure  at  a  value  of  0.1  R  where  R  = 
(a+b)/2,  and  again,  R  and  the  crack  length  are  held  at  a  value  of  one.  For  a  circular 
inclusion,  this  distance,  d,  is  0.1  times  the  radius  of  the  circle.  Note  that  for  a 
relatively  flat  ellipse,  or  for  a  value  of  m  greater  than  0.5,  the  location  of  the  close  tip 
of  the  crack  varies  rapidly  with  small  changes  in  the  orientation  angle  a.  This  effect 
is  shown  schematically  on  the  Figure.  In  effect,  once  the  crack  is  rotated  slightly  from 
the  x-axis,  it  is  shielded  by  the  flat  side  of  the  ellipse.  This  phenomenon  is  most 
striking  at  the  higher  ellipticity  ratios  (m=.9). 

Figure  (3.6)  presents  an  interesting  result  of  the  shielding  effect  of  a  flat 
ellipse.  In  this  Figure,  the  ellipticity  ratio  is  held  constant  at  a  value  of  0.9,  and  the 
distance  of  the  close  tip  of  the  crack  from  the  inclusion  is  varied  for  three  different 
angles.  The  angles  shown  are  very  small,  not  exceeding  4  degrees,  showing  that  very 
small  changes  in  orientation  of  the  crack  can  have  very  drastic  changes  in  the  stress 
Field  at  the  close  tip  of  the  crack.  There  is  a  30  percent  reduction  in  the  stress 
intensity  factor  for  an  angular  rotation  of  only  2  degrees,  if  the  crack  is  shielded  by 
the  ellipse.  As  the  close  tip  of  the  crack  is  moved  away  from  the  inclusion,  it  becomes 
unshielded,  and  asymptotically  approaches  the  stress  intensity  factor  for  a  crack  at  an 
orientation  angle  of  0  degrees  with  respect  to  the  inclusion  (see  Figure  3.2).  This 
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ellipticical  inclusions  (load  normal  to  crack,  R=l.,  L/R=l.,  k=1.67). 


shielding  is  an  expected  result,  and  should  provide  the  impetus  for  deflection  of  the 
path  of  a  naturally  growing  crack.  That  is,  of  course,  the  final  result  of  this  research, 
and  it  will  be  discussed  in  that  section  of  this  document. 

In  Figure  (3.7),  a  different  approach  was  taken  with  regard  to  the  location 
of  the  tip  of  the  crack  which  is  closest  to  the  inclusion.  In  this  figure,  the  distance 
held  fixed  is  that  from  the  close  tip  of  the  crack  to  the  origin  rather  than  to  the 
elliptical  inclusion.  Again,  the  crack  is  oriented  radially  with  respect  to  the  origin,  and 
normalized  stress  intensity  factors  are  plotted  for  three  different  ellipticity  ratios  versus 
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ratio  of  0.9  and  small  alpha  (load  normal  to  crack,  R=l.,  L/R=l.,  k=1.H). 


the  angle  that  the  crack  makes  with  the  x-axis.  The  distance  from  the  close  tip  of  the 
crack  to  the  origin  is  held  at  a  value  of  2.1R,  so  that  the  elliptical  inclusion  with  an 
m  of  1.0  (line  inclusion),  with  the  crack  at  an  angle  of  0  degrees  will  give  a  distance 
from  the  crack  to  the  inclusion  of  0.1R.  Note  that  there  is  still  shielding  of  the  crack, 
but  since  the  close  tip  of  the  crack  is  at  a  much  greater  distance  from  the  ellipse  than 
in  the  previous  figures,  the  shielding  effect  is  much  less  pronounced. 

Figure  (3.8)  is  the  only  one  in  which  the  crack  is  not  oriented  radially  with 
respect  to  the  origin.  In  this  figure  the  crack  is  either  horizontal  (Theta  =  0.)  or 
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Figure  3.7.  Stress  intensity  factor  versus  alpha  for  varying  ellipticical  inclusions 
with  crack  constant  distance  from  origin  (load  normal  to  crack,  R=l.,  L/R=l., 
k=1.67). 


vertical  (Theta  =  90.),  and  the  normalized  stress  intensity  factor  is  plotted  against  the 
x-coordinate  of  the  close  tip  of  the  crack.  The  y-coordinate  of  the  close  tip  is  held  at 
a  value  of  0.5R,  and  the  ellipse  has  an  ellipticity  ratio  of  0.8.  A  similar  result  was 
presented  for  the  circular  inclusion  by  Erdogan,  Gupta,  and  Ratwani  [7]  for  a  circular 
inclusion,  and  if  the  ellipticity  ratio  is  held  at  a  value  of  0.0,  the  present  technique 
matches  their  results.  Note  in  this  figure  that  the  crack  becomes  unshielded  at  an  x- 
coordinate  of  approximately  1.8R,  as  one  would  expect,  as  that  location  is  the  right 
edge  of  the  elliptical  inclusion. 
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horizontal  cracks  (theta  =  0  or  90,  load  normal  to  crack,  R=l.,  L/R=l.,  k=1.67). 


The  last  figure  of  the  series.  Figure  (3.9)  shows  the  effect  of  varying  the 
direction  of  the  applied  load  with  respect  to  the  crack.  In  this  figure,  the  distance 
between  the  inclusion  and  the  crack  was  held  fixed  at  d=.lR.  The  orientation  and 
elevation  of  the  crack  was  also  held  fixed,  at  0=0,  and  the  crack  aligned  along  the  x- 
axis.  The  direction  of  the  applied  load  on  the  crack  was  varied  from  a  value  of  0 
degrees  to  90  degrees,  for  two  ellipticity  ratios,  m=.4  and  .7,  and  for  two  different 
inclusions,  a  rigid  inclusion,  where  K  is  greater  than  1,  and  an  open  hole,  where  K=-l. 
Note  that  the  results  for  interaction  with  an  open  hole  are  significantly  different  than 
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those  for  a  rigid  inclusion,  but  that  the  t*  end  with  changing  load  angle  is  the  same  for 
all  curves.  The  Mode  I  stress  intensity  factor  increases  to  a  maximum  when  the 
applied  load  is  normal  to  the  crack,  as  would  be  expected.  This  figure  also  portrays 
another  interesting  phenomenon,  and  a  limitation  to  the  numerical  treatment  of  the 
problem.  As  the  angle  between  the  load  and  the  crack  approaches  zero,  note  that  the 
Mode  I  stress  intensity  factor  does  not  identically  approach  zero  for  three  of  the  four 
cases  shown  in  the  figure.  In  fact,  for  two  of  the  cases,  the  stress  intensity  factor 
becomes  less  than  zero.  A  negative  Mode  I  s^ess  intensity  factor  denotes 
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compression,  and  physically  represents  crack  face  overlap.  Mathematically,  what  is 
presented  is  correct,  but  physically,  the  two  crack  faces  cannot  overlap.  To  exclude 
overlap  in  the  numerical  procedure  would  necessitate  an  iterative  solution  procedure 
which  would  look  for  the  possibility  of  the  two  crack  faces  touching,  and  elimination 
of  that  portion  of  the  crack  from  the  solution.  That  procedure  has  not  been 
implemented  in  the  present  solution.  Also  note  that  there  is  one  case  in  which  the 
Mode  I  stress  intensity  factor  is  greater  than  zero  for  a  load  parallel  to  the  crack.  This 
is  both  physically  and  mathematically  correct,  as  the  presence  of  the  inclusion  can 
cause  a  tensile  stress  to  act  on  the  crack  even  when  the  load  is  in  the  same  direction 
as  the  crack. 

The  resubs  presented  in  this  section,  which  are  original,  both  verify  the 
solution  for  a  straight  crack  interacting  with  a  rigid  elliptical  inclusion,  and  also 
provide  some  interesting  motivation  for  the  prediction  of  the  path  of  the  crack  near  the 
inclusion.  These  results  are  valid,  however,  only  for  a  crack  which  is  straight.  They 
cannot  be  applied  to  a  crack  with  curvature,  because  as  the  crack  curves,  that  curvature 
provides  another  component  of  stress  which  is  not  accounted  for  in  the  present  method. 
In  essence,  the  problem  is  a  non-linear  one  in  that  the  crack  path  is  not  straight.  The 
method  presented  so  far  in  this  work  is  linear,  and  is  to  be  used  to  solve  a  non-linear 
problem.  The  next  chapters  in  this  work  will  describe  the  means  that  were  created  to 
allow  the  incorporation  of  this  non-linearity  into  the  present  solution,  and  therefore  to 
incrementally  predict  the  path  of  the  quasi-statically  growing  crack.  The  discussion 
and  the  results  that  are  presented  are  original. 
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Chapter  4 


CRACK  CURVATURE  AND  EXTENSION 


4.1  Crack  Kinking  Solutions 

As  described  in  the  survey  of  pertinent  literature,  there  are  several  solutions 
to  the  kinked  crack  problem.  Probably  the  one  most  referenced  recently,  and  most 
commonly  regarded  as  correct  is  the  one  published  by  K.  K.  Lo  [24]  in  1978.  His 
solution  unified  the  kinked  crack  problem  by  providing  correct  results  for  both  the 
finite  crack  extension,  and  an  infinitesimal  crack  extension.  At  the  time,  there  was  a 
good  deal  of  discussion  in  the  literature  about  this  problem,  as  described  in  the 
literature  review  portion  of  this  work,  and  there  was  also  a  good  deal  of  controversy. 
Since  that  time,  Lo’s  solution  has  been  accepted  as  being  correct,  and  has  been 
referenced  by  many  authors.  His  results  match  those  of  previous  authors  such  as 
Kitigawa  et  al.  [23]  for  the  finite  crack  extension,  and  Bilby  and  Cardew  [18]  for  the 
infinitesimal  crack  extension. 

The  solution  proposed  by  Lo  is  one  which  maps  the  kinked  crack  onto  a 
circle  using  a  conformal  mapping  technique.  He  then  calculates  potentials  for  the 
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crack  in  the  mapped  space,  solves  the  problem  in  the  mapped  space,  and  remaps  the 
solution  onto  the  kinked  crack.  His  mapping  is  conformal,  but  is  not  the  same 
mapping  as  that  used  in  the  present  solution  that  maps  an  elliptical  inhomogeneity  onto 
a  circular  region.  The  solution  proposed  by  Lo,  therefore,  cannot  be  used  directly  in 
the  present  work.  Lo’s  results,  however,  can  be  used  to  verification  the  results  of  the 
present  solution  once  that  solution  is  found. 

Another  interesting  solution  to  the  kinked  and  curved  crack  problem  is  that 
proposed  by  Sur  and  Altiero  [34],  In  their  solution,  as  described  in  the  review  of 
pertinent  literature  in  the  introduction,  the  unknown  function  is  the  crack  displacement 
itself  rather  than  the  derivative  of  that  displacement,  or  Burger’s  vector.  In  this  way, 
the  singularities  at  the  tips  of  the  crack  and  at  the  kink  are  avoided.  Results  are 
presented  for  this  method  which  are  in  excellent  agreement  with  other  published 
results. 


4.2  First  Order  Perturbation  Solution 

In  the  early  1970’s,  a  group  of  Russian  mathematicians  [25,  26,  and  27] 
solved  the  problem  of  a  crack  with  a  slightly  curved  extension,  as  shown  in  Figure  4. 1 . 
The  straight  portion  of  their  crack  extended  to  the  right  of  the  origin,  and  the  curved 
extension  extended  to  the  left,  as  shown  in  the  figure.  In  a  sense,  this  solution  avoids 
the  problems  of  the  weak  singularity  at  the  kink  by  assuming  that  the  crack  is 
smoothly  curved,  and  therefore  is  not  singular  at  the  kink.  They  used  the  complex 
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Figure  4.1  Schematic  of  crack  problem  solved  by  First  Order  Perturbation. 


potential  methods  of  Muskhelishvili,  and  wrote  an  integral  equation  that  describes  the 
singular  and  non-singular  parts  of  their  solution.  Their  solution  requires  knowledge 
of  the  local  spatial  derivatives  of  the  crack,  as  a  result  of  using  perturbation  methods. 
What  results  is  an  integral  equation  in  exactly  the  same  form  as  that  used  in  the 
present  method,  with  the  singular  and  non-singular  parts  separated  from  each  other, 
and  with  a  forcing  function  applied  along  the  crack  on  the  right  hand  side  of  the 
equation.  Their  integral  equation  can  be  restated  as  follows: 


J  —  T^*/t  +  J  R(s,x)b(x)dx  =  o*(s) 

-i S-T  -i 


(4.1) 


where  the  integral  is  given  on  the  interval  from  - 1  to  1  as  is  done  in  most  numerical 
crack  solution  methods,  including  Gerasoulis’  method,  and  the  star  on  the  forcing 
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function  denotes  that  the  stress  is  normalized  and  also  given  on  the  interval  from  -1 
to  1.  The  elements  of  their  kernel,  denoted  as  R(s,x)  arc  given  as  follows: 


R  _  Y\y*2  -  1  -  y7(*)y*(3+y*2)]  _  j 
"  "  (1  +  y*2)2 

R  =  y(y*2-i)(y*-y;C»» 

*  (1  +  y*2)2 

»  =  Yly'Q+y*2)  *  y\x)(y'2-i)] 

**  (i  +  y*2)2 

=  r(y*2-l)(l +/(*»*)  +  j 
^  (1  -  y*2)2 

where  the  capital  letters  and  stars  have  significance  as  follows: 


(4.2) 


>■•(*, T)  = 

s-x 

ns,T)  = 

s-x 


(4.3) 


in  that  they  are  approximations  to  the  first  and  second  derivatives  of  the  y-coordinate 
along  the  crack.  The  prime  notation  is  understood  to  be  differentiation  with  respect 
to  the  crack.  The  n  and  t  subscripts  refer  to  the  local  normal  and  tangential  directions 
along  the  crack,  and  the  x  and  y  subscripts  refer  to  the  x  and  y  components  of  the 
Burger’s  vector.  Close  inspection  of  the  four  elements  of  this  kernel  reveal  that  the 
quantities  in  the  numerator  of  each  term  are  functions  of  the  spatial  derivatives  and 
approximations  of  the  spatial  derivatives  taken  along  the  length  of  the  crack,  and  the 
denominator  is  the  square  of  the  length  of  the  curved  crack.  What  this  set  of  equations 


46 


contains,  therefore,  is  an  additional  set  of  terms  which  take  into  account  the  curvature 


of  the  crack,  and  can  be  directly  superposed  upon  any  straight  crack  solution  that  has 
a  similar  form,  to  mimic  the  effect  of  the  weak  singularities  inherent  in  crack  kinking 
and  curving.  In  the  present  work,  the  stresses  that  can  be  calculated  with  this  set  of 
equations  are  merely  superposed  upon  the  straight  crack  solution  for  a  curved  crack. 

However,  before  the  method  can  be  applied  to  the  current  solution  with  an 
elliptical  inhomogeneity,  local  spatial  first  and  second  derivatives  along  the  crack  must 
be  calculated  for  use  in  Equation  4.2. 


4.3  Parameterization  of  Crack  -  Cubic  Spline 

The  simplest  method  of  providing  these  derivatives  is  to  borrow  a 
technique  from  the  design  of  curves  and  surfaces,  the  notion  of  a  spline.  A  spline  is 
merely  a  set  of  smoothed  polynomial  interpolation  functions  which  are  used  to  evaluate 
a  function  at  any  point  along  an  interval,  if  the  function  is  only  known  at  discrete 
points  on  that  interval.  The  functions  used  are  ..nooih  in  the  sense  that  the  first 
derivative  is  continuous  along  the  interval,  and  are  independent  of  where  the  function 
is  known  or  is  to  be  evaluated.  Splines  have  been  used  intensively  in  Computer  Aided 
Design  for  the  description  of  space  curves  and  surfaces.  A  quadratic  spline  uses  a 
second  order  polynomial  to  describe  a  curved  line  segment,  a  cubic  spline  uses  a  third 
order  polynomial,  and  other  splines  use  other  orders  and  types  of  polynomials.  By  far 
the  simplest  splines  are  quadratic  and  cubic  splines.  Since  a  spatial  second  derivative 
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along  the  crack  is  needed,  the  cubic  spline  is  the  simplest  choice,  as  the  cubic 
polynomial  which  describes  the  line  segment  has  four  independent  constants  which  can 
be  varied.  This  is  the  correct  number  for  the  boundary  conditions  that  are  imposed  on 
the  endpoints  of  the  spline.  These  conditions  are  that  each  endpoint  of  an  interval 
must  be  connected  to  the  endpoint  of  the  next  interval,  and  the  first  derivative  must 
be  continuous  from  one  interval  to  the  next,  or  there  must  be  a  smooth  transition  from 
one  interval  to  the  next  with  no  kinks  in  the  graph  of  the  function.  This  lack  of  kinks 
in  the  graph  of  the  function  makes  the  spline  a  good  match  for  the  first  order 
perturbation  solution  which  also  avoids  kinks.  Therefore  there  are  two  conditions  at 
each  end  of  each  interval,  producing  a  total  of  four  conditions  which  must  be  satisfied. 
There  are  other  types  of  splines  which  are  relatively  difficult  to  program,  but  that  can 
describe  extremely  complex  shapes.  These  spline  functions  are  not  necessary  for  this 
work,  as  the  crack  path  will  be  reasonably  smooth.  Fortunately,  there  are  subroutines 
that  are  widely  available  for  fitting  a  cubic  spline  to  a  set  of  x  and  y  points.  The 
spline  routine  in  [44]  is  a  very  good  example.  The  algorithm  which  parameterizes  the 
crack  for  the  solution  provides  a  set  of  x  and  y  coordinates  at  the  integration  and 
collocation  points  along  the  crack.  These  locations  are  input  directly  into  the  spline 
routine,  which  calculates  the  second  derivative  at  each  location  where  there  is  a  point. 
The  distribution  of  first  derivatives  is  calculated  by  assessing  the  local  tangent  to  the 
crack  at  each  of  the  integration  and  collocation  points.  These  values  are  then  sent  to 
a  subroutine  which  uses  them  to  calculate  the  terms  in  equations  4.2. 
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4.4  Crack  Extension 


Once  these  results  have  been  generated,  the  algorithm  for  extending  the 
crack  is  relatively  simple.  At  each  integration  point  and  collocation  point,  there  is  a 
local  value  of  z,  or  a  local  position  along  the  crack,  which  has  x  and  y  components. 
There  is  also  a  local  value  of  the  angle  of  the  crack  4>,  which  is  initially,  when  the 
crack  is  straight,  equal  to  0,  the  initial  crack  angle.  There  is  also  a  set  of  first  and 
second  derivatives  along  the  crack  which  are  calculated  using  the  spline  routine  as 
described  above.  There  is  also  a  parameter  n  which  describes  the  size  of  the  matrix, 
and  the  number  of  integration  and  collocation  points.  There  are  2n  integration  points 
and  2n+l  collocation  points  which  describe  the  crack.  Since  the  computer  will  not 
allow  integers  to  be  incremented  by  a  value  of  0.5,  at  some  increment  of  crack 
extension,  n  must  increase  by  1,  as  this  is  the  basis  for  describing  the  size  of  the 
problem.  This  really  is  a  consequence  of  using  the  method  of  Gerasoulis,  as  his 
equations  are  based  upon  this  parameter,  n.  Incrementing  n  by  1  will  add  two 
integration  points  and  two  collocation  points  to  the  crack  parameterization.  What  is 
needed,  then,  is  to  decide  how  to  add  these  points  to  the  crack  parameterization.  The 
locations  of  the  integration  and  collocation  points  are  stored  in  complex  arithmetic 
arrays  in  the  computer,  in  order  of  where  they  are  located  along  the  crack.  The  values 
at  the  beginning  of  the  arrays  are  at  the  close  tip  of  the  crack,  and  the  values  at  the 
ends  of  the  arrays  are  at  the  far  tip  of  the  crack.  Since  it  is  growth  at  the  close  tip  of 
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the  crack  that  is  of  interest,  all  of  the  values  of  the  crack  parameterization  (integration 
and  collocation  points)  need  to  be  shifted  down  two  locations  in  their  respective  arrays 
on  the  computer.  This  is  done  by  starting  at  the  end  of  the  arrays,  adding  2  to  the 
maximum  extent  of  the  array,  and  copying  the  value  from  two  locations  before  that 
location  into  the  current  location.  In  this  way,  all  of  the  old  values  are  shifted  down 
2  spots  in  the  arrays.  Two  new  values  can  then  be  added  at  the  beginning  of  the 
arrays  for  the  z  location  of  the  new  close  crack  tip,  for  the  local  angle  0,  and  for  the 
local  first  and  second  derivatives.  This  increment  of  crack  extension  will  therefore  be 
of  length  1/n,  and  can  be  at  any  reasonable  angle  chosen.  This  technique  is  similar  to 
incrementing  pointers  in  the  Pascal  and  C  computer  languages.  Fortran,  unfortunately, 
does  not  share  this  data  type  with  these  two  other  languages,  and  slightly  more  coding 
is  required  to  implement  this  array  shifting. 


4.5  Results  and  Comparison  to  Previous  Work 

The  results  of  the  implementation  of  the  solution  described  above  are  given 
in  Figures  4.2  and  4.3.  These  figures  give  results  respectively  for  the  Mode  I  and 
Mode  II  stress  intensity  factors,  in  an  infinite  medium  with  no  inclusion,  for  an 
asymmetrically  kinked  crack  for  three  different  angles.  The  crack  is  asymmetrically 
kinked  in  that  it  has  a  kink  at  only  one  end,  as  shown  on  the  figures.  These  results 
compare  favorably  to  Kitigawa  et  al.  [23]  and  to  the  more  recent  results  of  Sur  and 
Altiero  [34].  The  results  given  by  Kitigawa  et  al.  were  verified  by  Lo  [24],  and  were 
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Figure  4.2.  Normalized  Mode  I  stress  intensity  factor  versus  normalized  kink  length 
for  three  different  kink  angles,  with  no  inclusion. 


used  as  a  verification  of  Sur  and  Altiero’s  work.  Note  that  the  results  of  the  present 
study  are  in  good  agreement  with  those  previously  published  results,  especially  for 
small  kink  angles.  This  is  a  very  encouraging  result,  as  the  formulation  using  the  first 
order  perturbation  solution  and  the  cubic  spline  crack  description  does  not  explicitly 
handle  the  singularity  at  the  kink.  In  fact,  the  present  solution  is  really  for  a  curved 
crack  rather  than  for  a  kinked  crack.  One  can  see,  however,  that  the  results  of  the 
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Figure  4.3.  Normalized  Mode  II  stress  intensity  factors  versus  normalized  kink 
length  for  three  different  kink  angles,  with  no  inclusion. 


present  solution  match  the  published  kinked  crack  solutions  quite  well  for  kink  angles 
less  than  about  30  degrees.  The  results  for  a  kink  angle  of  45  degrees  are  not  in  as 
good  agreement,  and  in  fact  fail  to  capture  the  singularity  as  the  kink  length  becomes 
small  with  respect  to  the  crack  length.  Fortunately,  it  is  expected  that  the  curvature 
of  a  naturally  growing  crack  will  be  smooth,  and  therefore  that  any  kink  angles  will 
be  small.  The  first  order  perturbation  solution,  therefore,  seems  to  give  good  results. 
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and  allows  the  approximation  of  the  weak  singularities  due  to  crack  kinking  without 
explicitly  describing  that  weak  singularity. 
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Chapter  5 

CRACK  PATH  PREDICTION 


5.1  Review  of  Crack  Propagation  Laws  and  Selection  of  Proper  Criterion 

There  are  several  means  which  have  been  proposed  to  predict  the  direction 
of  the  onset  of  crack  propagation.  Most  of  these  fall  into  two  categories,  stress 
methods  and  energy  methods.  Several  authors  [36,  37,  38,  and  17]  have  proposed  a 
criterion  for  growth  which  is  called  the  Maximum  Normal  Stress  Criterion,  in  which 
one  looks  for  the  direction  where  the  Mode  I  stress  intensity  factor  is  greatest,  and 
concludes  that  the  crack  propagates  in  that  direction.  Other  authors  [39,  40,  41,  and 
43]  have  proposed  what  has  come  to  be  called  the  Maximum  Strain  Energy  Release 
Rate  Criterion,  in  which  one  calculates  the  direction  at  which  the  maximum  energy  is 
released,  and  then  concludes  that  the  crack  propagates  in  that  direction.  A  very  good 
review  of  those  two  major  classes  of  crack  propagation  criteria,  and  some  less  well 
known  ideas,  is  given  by  Palinaswamy  and  Knauss  [17].  These  authors  correctly  state 
that,  for  isotropic  materials,  and  for  cracks  whose  kink  angle  is  fairly  small  (less  than 
about  30  degrees),  the  two  criteria  produce  identical  results.  Actually,  there  has  been 
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a  significant  amount  of  controversy  over  which  criterion  is  correct.  The  two  criteria 
begin  to  diverge  when  the  angle  between  the  load  and  the  crack  is  greater  than  30 
degrees.  There  have  been  several  experiments  performed  in  attempts  to  settle  this 
controversy,  with  limited  results.  The  divergence  between  the  two  criteria  is  a 
maximum  when  the  load  is  parallel  to  the  crack,  a  condition  where  the  Maximum 
Stress  theory  predicts  a  kink  angle  of  70.5  degrees,  whereas  the  Maximum  Strain 
Energy  Release  Rate  criterion  predicts  an  angle  of  77.4  degrees.  Several  authors  have 
produced  figures  which  portray  this  divergence,  and  also  show  the  scatter  of 
experimental  data  (Williams  and  Ewing  [39],  Finnie  and  Saith  [40],  Ewing  and 
Williams  [41],  and  Wu  [42]).  The  controversy  has  not  yet  been  settled.  For  this  work, 
however,  the  fact  that  it  is  expected  that  the  curvature  of  a  naturally  growing  crack 
will  not  exhibit  even  a  30  degree  kink  leads  to  the  conclusion  that  either  criterion  will 
be  acceptable  for  the  present  discussion. 

What  results  then,  is  the  selection  of  which  criterion  best  suits  the  crack 
path  prediction  algorithm.  The  strain  energy  released  by  an  increment  of  crack  growth 
can  be  calculated  as  follows: 
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where  both  the  energy  and  the  stress  intensity  factors  are  normalized  to  the  value 
corresponding  to  the  stress  intensity  of  a  finite  crack  in  an  infinite  medium  that  is  not 
under  the  influence  of  an  inclusion.  Since,  in  the  direction  of  crack  growth,  the  value 
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of  the  Mode  II  stress  intensity  factor  is  nearly  zero,  and  the  value  of  the  Mode  I  stress 
intensity  factor  changes  very  little,  there  will  be  very  small  changes  in  the  strain 
energy  calculated  as  the  square  of  a  number  close  to  one  added  to  the  square  of  a 
number  close  to  zero.  In  fact,  for  a  change  in  direction  of  as  much  as  5  degrees,  the 
strain  energy  changes  less  than  1  percent.  On  the  computer,  looking  for  a  local 
maximum  of  a  number  that  is  changing  by  about  0.1  percent  per  iteration  is  a  very 
difficult  task. 

The  application  of  the  Maximum  Normal  Stress  Criterion,  requires  that  the 
Mode  II  stress  intensity  factor  be  made  to  come  as  close  to  zero  as  possible.  This 
stress  intensity  factor  may,  in  some  situations,  not  identically  equal  zero,  as  there  may 
be  a  Mode  II  component  that  is  present  in  any  case.  That,  in  fact,  is  the  situation  that 
arises  in  this  problem,  there  is  a  small  amount  of  Mode  II  loading  present  just  because 
of  the  presence  of  the  inclusion.  The  application  of  that  criterion,  therefore,  is 
implemented  as  a  search  for  the  direction  in  which  the  Mode  II  stress  intensity  factor 
is  an  absolute  minimum  in  magnitude.  It  is  reasonably  simple  to  find  an  absolute 
minimum  on  the  computer,  rather  than  a  relative  maximum.  For  this  reason,  really  one 
of  simplicity,  the  algorithm  which  predicts  the  pain  of  the  crack  searches  for  the  angle 
at  which  the  Mode  II  stress  intensity  factor  is  a  minimum  in  magnitude. 
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5.2  Path  Prediction  Algorithm 

It  would  appear,  using  the  Maximum  Normal  Stress  Criterion,  that  the  path 
prediction  would  be  a  reasonably  simple  task,  if  one  is  merely  to  find  the  direction  at 
which  the  shear  stress  vanishes.  Erdogan  and  Sih  [36]  in  fact,  provide  a  formula 
which  can  be  used  to  predict  this  initial  tip-off  angle  given  the  stress  intensity  factors 
at  the  original  crack  tip.  They  use  this  equation  to  predict  the  initial  tip-off  angle  of 
a  crack  which  is  loaded  in  pure  shear,  with  very  good  results.  This  procedure  was 
attempted  in  some  initial  tries  at  obtaining  a  path  prediction,  with  no  success.  In  the 
present  solution,  the  effect  of  the  elliptical  inclusion  is  to  alter  the  stress  field 
significantly  with  advance  of  the  crack.  Therefore,  the  angle  that  would  be  predicted 
from  the  initial  crack  stress  intensity  factors  is  not  the  angle  at  which  the  natural  crack 
will  grow.  With  each  infinitesimal  crack  extension  that  occurs  with  the  naturally 
growing  crack,  the  stress  field  changes  significantly  due  to  the  interaction  with  the 
inclusion.  This  phenomenon  is  illustrated  in  Figures  5.1  through  5.7.  Each  of  these 
figures  displays  contours  of  the  maximum  shear  stress  near  the  tip  of  a  growing  crack 
as  it  is  interacting  with  an  inclusion.  These  are  the  same  stress  contours  that  would 
appear  as  isoclines  in  a  photoelastic  stress  analysis  of  this  problem.  The  ellipticity 
ratio  of  the  inclusion  is  0.7,  and  the  initial  close  crack  tip  is  located  at  a  z  value  of 
(2. 5, .3),  or  2.5  in  the  x-direction,  and  0.3  in  the  y-direction.  The  initial  value  of  n  is 
5,  so  the  increments  of  crack  growth  are  0.2  in  length.  The  ellipse  has  a  unit  R,  and 
the  crack  is  initially  1.0  units  in  length.  Note  that  the  crack  grows  until  it  reaches  the 
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growth  (load  normal  to  crack,  0=0.,  R=l.,  L/R=l.,  k=1.67). 

ellipse,  and  then  is  stopped.  This  is  an  actual  prediction  of  crack  path,  and  will  be  one 
that  is  presented  in  a  figure  later  in  this  work.  Careful  investigation  of  these  seven 


figures  shows  that  the  contours  of  stress  tilt  significantly  with  the  advance  of  the  crack. 
Also,  they  are  not  even  similar  in  form  on  either  side  of  the  crack.  In  the  solution 


presented  by  Erdogan  and  Sih  [36],  these  contours  of  stress  were  identically  reflected 
on  opposite  sides  of  the  crack.  Therefore,  instead  of  a  linear  crack  advance  problem, 
the  present  situation  is  in  fact  non-linear.  The  significance  of  this  and  of  these  figures 
is  that  the  prediction  algorithm  that  is  developed  can  at  most  be  incrementally  linear. 
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and  the  correct  crack  propagation  direction  cannot  be  determined  with  information 
available  before  the  crack  actually  propagates. 

The  prediction  algorithm  for  this  problem  therefore  must  be  iterative  in 
nature.  An  initial  guess  is  made  at  the  angle  at  which  the  crack  is  going  to  grow,  and 
then  is  checked  to  see  if  that  guess  is  correct.  Actually,  the  prediction  algorithm  that 
is  used  to  predict  the  direction  of  crack  growth  is  reasonably  straight-forward.  It 
begins  by  initially  solving  the  problem  of  a  straight  crack  with  unit  length,  as 
described  previously.  Then  an  increment  of  crack  length  is  added  to  the  crack  by 
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incrementing  the  problem  size  parameter  n,  then  the  arrays  are  incremented  by  2,  and 
the  values  are  moved  down  2  locations.  Next,  a  choice  for  the  angle  at  which  the 
crack  will  deflect  from  its  present  path  is  made.  A  very  efficient  choice  is  for  the 
angle  to  grow  in  the  same  direction  that  the  crack  is  already  growing,  or  for  each 
increment,  for  there  to  be  no  more  deflection.  For  a  naturally  growing  crack,  we 
expect  that  the  path  will  be  relatively  smooth,  or  that  there  will  be  very  little  deflection 
observed  during  each  increment  of  growth.  Therefore,  the  angle  that  is  implemented 
is  no  deflection  from  the  current  path,  or  algorithmically,  the  angle  calculated  during 


60 


the  last  increment  of  growth.  The  problem  is  then  solved  again,  and  stress  intensity 
factors  are  calculated.  Once  these  values  are  known,  the  program  proceeds  to  the  next 
iteration.  At  this  iteration,  a  known  angle  is  added  to  the  initial  guess,  and  the 
problem  is  solved  again.  At  this  and  subsequent  iterations,  however,  there  is  no  need 
to  extend  the  crack,  nor  to  increment  the  arrays  nor  move  the  values  in  the  arrays. 
The  first  iteration,  at  zero  angle  and  one  increment  of  crack  growth  has  done  the 
extension  of  the  crack,  and  the  rearranging  of  the  values  in  the  arrays.  To  perform  the 
second  and  subsequent  iterations  for  this  increment  of  crack  growth,  therefore,  the  first 
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Figure  5.5.  Contours  of  maximum  shear  stress  after  fourth  increment  of  crack 
growth  (load  normal  to  crack,  0=0.,  R=l.,  original  L/R=l.,  k=1.67). 


two  values  in  each  array  need  to  be  replaced  with  new  values,  as  the  new  crack  tip  is 
at  a  new  value  of  z,  with  a  new  angle,  and  new  first  and  second  derivatives.  The 
problem  then  continues  to  be  solved  again  and  again  with  new  kinked  tips.  The  angle 
chosen  for  the  initial  tip-off  after  the  first  iteration  with  no  deflection  is  an  angle  of 
4  degrees.  This  angle  was  chosen  after  quite  a  bit  of  experimentation  in  reducing  the 
number  of  iterations  required  to  converge  on  the  answer. 

The  next  step  is  to  determine  whether  the  new  increment  of  crack  growth 
is  in  the  correct  direction.  To  do  this,  the  new  Mode  II  stress  intensity  factor  that  is 
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calculated  in  the  current  iteration  is  compared  to  the  result  from  the  previous 
calculation.  If  the  new  Kn  is  less  than  the  old  value,  the  crack  tip-off  angle  is 
incremented  by  another  4  degrees  and  the  problem  is  solved  again.  If  the  new  Mode 
II  stress  intensity  factor  is,  however,  greater  than  the  old  value,  the  algorithm  must 
have  overshot  the  correct  angle,  or  gone  the  wrong  direction.  Since  it  is  not  known 
for  certain  whether  there  has  been  an  overshot  or  the  prediction  has  gone  the  wrong 
direction,  the  best  thing  that  can  be  done  is  to  cut  the  increment  in  half  and  go  back 
the  other  direction.  If  there  were  human  interaction  within  this  algorithm,  the  human 
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could  guess  at  the  new  value  of  the  kink  angle,  and  tell  the  computer  to  try  that  value. 
The  intent  of  this  algorithm,  however,  was  to  make  the  prediction  automatic,  and  not 
require  human  intervention.  It  needs,  therfore,  to  converge  on  the  correct  angle  with 
a  minimum  of  iterations,  and  no  human  intervention.  The  new  guess  at  0  will 
therefore  be  2  degrees  less  than  the  old  value  of  the  crack  kink  angle.  Again,  once  the 
problem  is  solved,  there  is  a  new  value  for  Kn  with  which  to  compare  the  old  value. 
The  guessing  and  going  backwards  with  half  the  increment  continues  for  one  more 
reversal  in  direction,  until  the  increment  in  angle  is  1  degree.  Once  the  increment  is 
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reduced  to  1  degree,  and  a  minimum  has  been  found  again,  the  kink  angle  is  defined 
as  converged.  This  provides  a  solution  that  is  incorrect  by  at  most  half  that  amount. 
That  is,  the  correct  answer  will  be  off  by  at  most  1/2  of  a  degree  if  the  increment  of 
change  to  the  guess  is  1  degree. 

Once  the  algorithm  has  converged  on  an  acceptable  answer,  the  computer 
program  goes  back  with  its  final  answer  for  the  crack  kink  angle,  0,  and  solves  the 
new  problem  with  the  new  angle.  The  prediction  continues  with  the  new  crack  and 
the  procedure  begins  again  for  the  next  increment.  There  is  only  one  thing  that  must 
be  checked  at  each  increment  of  growth,  whether  or  not  the  crack  has  grown  into  the 
inclusion.  This  is  also  reasonably  simple.  The  distance  of  the  new  crack  tip  from  the 
origin  can  be  easily  calculated,  as  can  the  angle  that  a  vector  from  the  origin  to  that 
location  would  make  with  respect  to  the  x-axis.  The  distance  to  the  edge  of  the  ellipse 
at  that  angle  can  then  be  calculated,  and  a  comparison  made  of  the  two  values.  If  the 
edge  of  the  ellipse  is  closer  than  the  crack  tip,  the  crack  is  still  outside  the  inclusion, 
and  the  algorithm  can  proceed.  If  not,  the  crack  has  grown  to  the  edge  of  the  ellipse, 
and  the  program  ends  the  prediction.  Actually,  the  crack  will  grow  into  the  inclusion, 
or  along  the  edge  of  the  inclusion,  but  due  to  the  limitations  of  rigidity  of  the 
inclusion,  these  problems  cannot  be  solved  presently. 
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5.3  Verification  of  Prediction 


Sumi  et  al.  [33]  performed  finite  element  crack  path  predictions  near  an 
open  hole,  and  also  performed  experiments  to  verify  their  predictions.  Figure  5.8 
shows  Sumi’s  predictions  and  experiments,  and  predictions  done  with  the  present 
algorithm.  The  inclusion  in  the  present  solution  is  a  rigid  one,  but  the  value  of  k 
(Muskhelishvili’s  constant)  can  be  set  to  any  value  required.  As  stated  previously, 
Muskhelishvili  noted,  in  his  work,  at  a  value  of  K  =  -1  the  solution  for  a  rigid 
inclusion  reduces  to  the  solution  for  a  hole.  This  is  helpful  in  that  it  allows  the 
calculation  of  the  interaction  of  a  crack  with  an  open  hole  to  be  done  using  the  results 
for  a  rigid  inclusion.  This,  in  fact,  is  what  was  done  to  produce  the  predictions  shown 
in  the  figure.  Note  that  the  predictions  done  with  the  current  algorithm  match  the 
experimental  results  reasonably  well,  especially  when  the  crack  passes  the  open  hole. 

Figure  5.9  depicts  the  path  of  a  crack  which  is  growing  naturally  near 
either  an  open  hole,  a  rigid  inclusion,  or  no  inclusion  at  all.  Note  that  the  hole  attracts 
the  crack,  as  has  been  observed  by  experiments  (Sumi  [33]),  the  rigid  inclusion 
deflects  the  path  of  the  crack  as  would  be  expected,  and  the  case  with  no  inclusion  has 
no  effect  on  the  crack  at  all,  as  would  be  expected.  All  three  of  these  cases  are  actual 
predictions  using  the  method  that  is  described  in  this  work. 

It  appears  that  the  predictions  of  a  crack  path  are  reasonable.  The  results 
using  the  present  method  compare  well  to  predictions  and  experiments  done  by  others. 
What  remains,  then,  is  to  exercise  the  solution  in  order  to  learn  about  the  paths  of 
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Figure  5.8.  Verification  of  prediction  by  comparison  to  experimental  results  (load 
normal  to  crack,  R=l.,  original  L/R=l.,  K=-l.). 


cracks  near  inclusions. 

5.4  Results 

Figures  5.10  through  5.14  depict  results  of  the  crack  path  prediction.  All 
of  these  figures  represent  the  actual  path  of  the  cracks.  In  each  figure,  the  initial  crack 
is  not  shown,  but  is  of  length  1.0/R,  extending  in  the  positive  x  direction  from  the 
right-most  end  of  the  path  shown  in  the  figure.  The  x  and  y  coordinates  shown  in 
each  of  the  figures  represent  distances  from  the  origin  normalized  by  the  average 
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Figure  5.9.  Paths  of  cracks  given  an  open  hole,  no  inclusion,  and  a  rigid  circular 
inclusion. 

radius  of  the  ellipse,  R.  There  is  therefore  no  labelling  in  any  of  these  figures  of  the 
axes,  as  they  are  intended  to  portray  x  and  y  physical  space.  Figure  5.10  is  the  result 
of  a  parametric  study  on  the  ellipticity  ratio  of  the  inclusion.  In  this  figure  the  initial 
location  of  the  close  tip  of  the  crack  is  at  an  x-coordinate  of  2.5,  and  a  y  coordinate 
equal  to  the  uppermost  extent  of  the  edge  of  the  inclusion.  For  instance,  for  an 
inclusion  with  an  ellipticity  ratio  of  0.6,  the  initial  close  tip  of  the  crack  is  at  x  =  2.5, 
and  y  =  .4.  Note  that  the  inclusions  with  ellipticity  ratios  greater  than  0.4  tend  to 

68 


Figure  5.10.  Paths  of  cracks  initially  aligned  with  the  top  edge  of  inclusions  with 
varying  ellipticity  ratios  (load  normal  to  crack,  R=l„  original  IVR=1.,  k=1.67). 

attract  the  crack,  and  the  crack  does  not  escape  the  inclusion.  For  an  ellipticity  ratio 
of  0.4  and  less,  the  crack  is  deflected,  and  escapes  the  inclusion.  In  this  figure,  the 
ellipse  is  rigid.  Figure  5.11  portrays  exactly  the  same  parametric  study,  only  in  this 
case,  the  value  of  k  is  set  to  -1,  modelling  the  interaction  of  the  cracks  with  open 
elliptical  holes  with  varying  ellipticity  ratios.  Note  in  this  figure  that  the  crack  never 
escapes  the  hole. 

Figure  5.12  portrays  a  slightly  different  parametric  study  in  which  the 
initial  location  of  the  crack  is  held  fixed,  at  x  =  2.5  and  y  =  .6,  as  is  the  ellipticity 
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normal  to  crack,  R=l.,  original  L/R=l.,  K--1.) 


ratio,  at  a  value  of  0.4,  and  the  initial  angle  that  the  crack  makes  with  the  x-axis  is 
varied  from  0  to  30  degrees  in  increments  of  5  degrees.  Note  that  the  crack  only 
escapes  being  attracted  into  the  side  of  the  inclusion  for  the  case  of  a  0  degree 
(horizontal)  crack.  This  is  the  only  case  in  which  a  line  drawn  from  the  tip  of  the 
crack,  aligned  parallel  to  the  initial  crack  would  also  not  interfere  with  the  inclusion. 

In  Figures  5.13  and  5.14,  the  ellipticity  ratio  is  again  held  fixed  for  each 
figure,  as  is  the  initial  angle  of  the  crack.  In  these  figures,  the  initial  y-coordinate  of 
the  close  tip  of  the  crack  is  varied,  to  study  the  extent  to  which  an  inclusion  will 
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Figure  5.12.  Paths  of  cracks  with  varying  initial  angles  with  respect  to  the  x-axis 
(load  normal  to  crack,  R=l.,  original  L/R=l.,  k=1.67). 

attract  or  repel  a  crack.  Figure  5.13  depicts  crack  paths  for  a  rigid  inclusion  with  an 
ellipticity  ratio  of  0.9  (a  long,  flat  inclusion).  Note  that  cracks  with  a  y-coordinate  of 
0.4  escape  the  side  of  the  inclusion,  whereas  cracks  with  original  y-coordinate  of  0.3 
and  below  are  attracted  into  the  side  of  the  inclusion. 
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Figure  5.13.  Paths  of  cracks  interacting  with  inclusion  with  ellipticity  ratio  0.9,  and 
starting  at  varying  initial  y-coordinates  (load  normal  to  crack,  R=l.,  original  L/R=l., 
k=1.67). 


In  Figure  5.14,  the  ellipticity  ratio  of  the  inclusion  is  0.4.  For  this 
inclusion,  the  crack  must  have  an  initial  y-coordinate  of  0.6  before  escaping  the 
inclusion.  This  location  is,  however,  parallel  to  the  upper  edge  of  the  inclusion.  For 
the  previous  case,  cracks  were  attracted  toward  the  inclusion  from  a  y-coordinate 
considerably  above  the  upper  extent  of  the  inclusion.  Long,  flat  inclusions  therefore 
have  a  greater  tendency  to  attract  cracks  than  do  shorter  rounder  inclusions.  It  is  also 
interesting  to  look  at  the  tip  of  each  crack  in  Figure  5.14  as  it  approaches  the 
inclusion.  Cracks  approaching  the  inclusion  near  its  major  axis  are  merely  attracted 
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Figure  5.14.  Paths  of  cracks  interacting  with  an  inclusion  with  ellipticity  ratio  0.4, 
and  varying  y-coordinates  (load  normal  to  crack,  R=l.,  original  L/R=l.,  k=1.67). 

There  are  many  more  of  these  types  of  studies  which  can  be  done  to  learn 
interesting  and  useful  things  about  the  paths  of  cracks  near  inclusions  and  holes. 
These  results  are  presented  therefore  both  to  provide  the  reader  a  sense  of  the 
flexibility  and  general  nature  of  the  solution  as  well  as  providing  some  general  insight 
into  the  quasi-static  behavior  of  cracks  near  inclusions  and  holes. 
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5.3  Energy  Absorption  Due  to  Crack  Deflection 


There  are  two  components  to  the  additional  energy  absorption  associated 
with  a  non-straight  crack  path  near  an  inclusion.  The  first  component  is  the  effect  of 
the  local  shielding  of  the  crack  due  to  the  inclusion.  The  ratio  of  the  local  stress 
intensity  factors  to  the  applied  values,  or  in  the  present  case,  the  normalized  stress 
intensity  factors,  quantifies  the  amount  of  local  shielding  due  to  the  presence  of  the 
inclusion.  In  the  case  of  an  open  hole,  this  ratio  is  greater  than  one,  and  therefore 
describes  the  amount  of  anti-shielding  due  to  the  presence  of  the  hole.  This  effect  is 
what  is  calculated  and  has  been  presented  in  this  work.  The  energy  absorbed  or 
energy  shielding  ratio  can  be  calculated  by  determining  the  strain  energy  release  rate, 
as  follows: 


G_  =  *?  ^ 

G-  (om/na)2 


(5.1) 


where  K,  and  K„  and  G  are  the  local  values,  and  they  are  normalized  by  the  values 
that  would  exist  when  there  is  no  inclusion.  For  the  case  where  there  is  no  inclusion, 
therefore,  this  equation  gives  a  value  of  1.  In  the  presence  of  an  inclusion,  if  the  crack 
is  shielded  by  the  inclusion,  the  local  K,  will  be  less  than  1,  and  therefore,  the 
normalized  G  will  also  be  less  than  1.  This  decrement  is  energy  can  be  thought  of  as 
the  interaction  toughness  ratio,  or  the  ratio  of  energy  release  rate  for  an  inclusion 
present  to  that  with  no  inclusion.  In  the  present  work,  this  ratio  can  be  calculated  at 
each  increment  of  crack  growth,  and  the  total  can  be  averaged  by  the  number  of 
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increments  of  crack  growth  that  have  occurred. 

The  second  component  of  the  energy  absorption  is  that  due  to  the  path 
itself.  The  interaction  toughness  ratio  described  above  inherently  assumes  a  rectilinear 
crack  path.  However,  if  the  crack  is  deflected  from  a  rectilinear  path,  as  it  will  be 
near  an  inclusion,  then  the  path  itself  also  represents  a  toughening  mechanism.  This 
crack  deflection  toughness  ratio  can  also  be  assessed  using  the  technique  described  in 
this  work,  because  the  final  location  of  the  advancing  crack  tip  is  known,  and  the 
distance  that  the  crack  would  have  travelled  in  its  original  rectilinear  direction  is  also 
known. 

Therefore,  to  assess  the  overall  energy  absorption  due  to  the  crack 
interacting  with  an  inclusion,  and  being  deflected  by  the  inclusion  during  its  quasi¬ 
static  crack  growth,  both  of  these  ratios  should  be  used.  The  total  energy  absorbed, 
or  total  toughness  ratio  would  then  be  the  product  of  these  two  ratios.  An 
implementation  of  that  idea  would  be  as  follows  for  an  initially  horizontal  crack:  note 
the  final  x-location  of  the  growing  tip  of  the  initially  horizontal  crack,  keeping  track 
of  how  many  increments  of  crack  growth  have  taken  place;  then  divide  the  total  x- 
distance  that  the  crack  has  travelled  by  the  distance  that  it  could  have  travelled  if  it 
had  not  curved  (increment  length  times  number  of  increments).  This  number  will  be 
less  than  one,  and  will  also  represent  an  energy  decrement,  or  an  amount  of  energy 
absorbed  by  path  deflection  only.  This  path  deflection  ratio  can  then  be  multiplied  by 
the  interaction  toughness  ratio  calculated  as  an  average  over  all  of  the  increments  of 
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crack  growth,  to  determine  the  overall  energy  absorption  or  toughness  ratio  for  the 
crack. 

Energy  absorption  for  several  different  cracks  has  been  calculated  using 
this  procedure,  for  a  both  a  rigid  inclusion  and  an  open  hole,  and  is  presented  in 
Tables  5.1  and  5.2.  The  crack  parameters  are  the  same  as  those  for  the  paths 
calculated  for  Figures  5.10  and  5.11,  where  the  initial  tip  of  the  crack  is  at  an  x- 
coordinate  of  2.5/R,  and  the  y-coordinate  is  at  the  uppermost  extent  of  the  edge  of  the 
inclusion.  Ellipticity  ratios  varying  from  .9  to  0.  ere  displayed  in  the  table.  The 
numbers  in  the  table  represent  the  various  toughness  ratios  described,  where  the  second 
column  is  the  overall  toughness  ratio,  and  the  third  and  fourth  columns  are  the  path 
deflection  or  shortening  ratio,  and  the  shielding  or  interaction  ratio  respectively. 

It  is  interesting  to  note  from  a  first  glance  at  Table  5.1  that  the  greatest 
energy  absorption  occurs  when  the  crack  is  attracted  into  the  side  of  the  inclusion. 
However,  upon  closer  inspection  of  the  values  in  the  table,  one  can  see  that  it  is 
important  to  compare  similar  cracks.  The  energies  reported  in  the  table  for  ellipticity 
ratios  greater  than  0.4  are  for  cracks  which  are  considerably  shorter  than  those  for 
ellipticity  ratios  0.4  and  less,  since  they  are  attracted  into  the  side  of  these  inclusions. 

For  those  cracks  which  pass  the  inclusion,  the  shortening  ratio  is  relatively 
constant.  This  is  a  reasonable  result,  referring  to  Figure  5.10.  In  this  figure,  it  appears 
that  the  paths  of  the  cracks  which  escape  the  ellipse  are  nearly  parallel  to  each  other 
for  the  set  of  conditions  described.  It  is  expected,  therefore,  that  their  crack  length 
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Table  5.1.  Toughness  ratios  for  cracks  initially  parallel  to  the  top  of  the  inclusion, 
for  differing  inclusion  ellipticity  ratios. 


Ellipticity 

Ratio 

Toughness 

Ratio 

Crack  Length 
Shortening 

Shielding  Ratio 

.9 

.713 

.795 

.897 

.8 

.665 

.829 

.801 

.7 

.659 

.853 

.773 

.6 

.662 

.886 

.748 

.5 

.672 

.896 

.750 

.4 

.937 

.962 

.975 

.3 

.944 

.961 

.982 

.2 

.952 

.962 

.989 

.1 

.958 

.962 

.997 

0. 

.962 

.962 

1.00 

shortening  ratios  would  all  be  similar.  For  those  cracks  which  are  attracted  into  the 
side  of  the  ellipse,  all  of  the  shortening  ratios  will  be  different,  as  the  cracks  are  all 
of  different  lengths.  If  the  toughness  ratios  of  the  cracks  interacting  with  inclusions 
with  ellipticities  of  0.5  and  greater  are  compared,  it  would  seem  that  more  energy  is 
absorbed  by  shorter,  fatter  inclusions.  In  these  cases,  the  shielding  ratio  dominates  the 
energy  absorption.  For  inclusions  with  ellipticities  less  than  0.5,  however,  exactly  the 
opposite  is  true.  If  this  result  is,  in  fact,  the  true  case  in  real  materials,  it  could  have 


Table  5.2  Toughness  ratios  for  cracks  initially  parallel  to  the  top  of  the  hole,  for 
different  hole  ellipticity  ratios 


Ellipticity 

Ratio 

Toughness 

Ratio 

Crack  Length 

Ratio 

Anti-Shielding 

Ratio 

.9 

3.45 

.819 

4.21 

.8 

3.21 

.809 

3.97 

.7 

2.48 

.838 

2.96 

.6 

2.36 

.880 

2.68 

.5 

2.15 

.865 

2.48 

.4 

1.89 

.879 

2.25 

.3 

1.69 

.892 

1.89 

.2 

1.72 

.894 

1.92 

.1 

1.73 

.905 

1.91 

0. 

1.55 

.905 

1.71 

a  strong  impact  on  the  design  of  the  shapes  of  reinforcements. 

Table  5.2  shows  a  somewhat  different  result,  in  that  the  overall  toughness 
ratios  are  greater  than  one,  as  these  are  results  for  cracks  growing  toward  an  open  hole. 
Since  all  of  these  cracks  are  eventually  stopped  by  the  open  hole,  a  comparison  of  all 
of  the  results  together  in  this  table  is  reasonable.  Note  that  the  crack  length  shortening 
ratio  or  crack  deflection  toughness  ratio  remains  relatively  constant,  whereas  the  anti¬ 
shielding  ratio,  or  crack  inclusion  interaction  ratio  varies  significantly,  and  provides 
a  majority  of  the  propensity  for  crack  growth.  Also,  long,  thin  holes  have  a 
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considerably  stronger  attractive  effect  on  the  crack  than  do  rounder,  more  circular 
holes.  This  result  could  also  have  a  strong  effect  on  toughening  mechanisms  and 
reinforcing  mechanisms  in  materials,  as  it  seems  to  suggest  that  if  the  reinforcement 
becomes  debonded  from  the  surrounding  material,  the  resulting  open  hole  will  attract 
a  crack,  and  cause  it  to  propagate. 
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Chapter  6 
CONCLUSIONS 


6.1  Summary  of  Findings 

In  this  work,  a  new  technique  has  been  presented  which  will  predict  the 
path  of  a  naturally  growing  crack  near  a  rigid  elliptical  inclusion,  where  the  stress  field 
can  be  modelled  as  two  dimensional.  The  technique  uses  a  boundary  integral  approach 
to  the  solution  of  the  stress  field  at  the  tips  of  the  crack  which  is  in  general  not 
straight,  and  is  interacting  with  a  rigid  elliptical  inclusion.  The  crack  path  is 
parameterized  as  a  cubic  spline,  and  the  results  of  both  a  Green’s  Function  solution  to 
the  interaction  of  a  dislocation  with  an  elliptical  inclusion,  and  a  first  order 
perturbation  solution  to  account  for  the  generally  curvilinear  nature  of  the  crack  have 
been  employed.  The  singular  nature  of  the  stresses  is  accounted  for  using  a  numerical 
technique  which  describes  the  distribution  of  dislocations  along  the  crack  as  a 
piecewise  quadratic  polynomial  to  transform  the  problem’s  resulting  integral  equations 
into  algebraic  equations  well  suited  to  a  matrix-type  solution.  Results  of  each  step  of 
the  analysis  have  been  verified  with  previously  published  results  that  are  commonly 
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accepted  as  correct,  and  with  experimental  results  of  a  crack  propagating  near  an  open 
circular  hole.  New  and  significant  results  are  also  presented  as  paths  of  cracks 
interacting  with  inclusions  of  differing  ellipticity  ratios,  and  at  different  orientations 
with  respect  to  the  initial  crack.  The  technique  is  much  less  computationally  intensive 
than  current  combination  analytical  and  finite  element  approaches  to  crack  path 
prediction.  All  of  the  predictions  presented  in  this  work  were  generated  on  an 
engineering  workstation  type  computer.  There  was  no  need  to  go  to  a  larger  computer 
for  the  predictions.  This  indicates  that  sensitivity  studies  using  the  technique  are  quite 
easily  accomplished,  and  provide  some  new  insight  into  the  mechanisms  of  fracture 
near  inclusions  and/or  reinforcements  in  materials.  The  angle  that  the  crack  has  with 
the  x-axis  is  not  limited,  nor  is  the  angle  that  the  applied  load  makes  with  respect  to 
the  initial  crack,  except  for  the  case  in  which  the  crack  faces  would  overlap.  Original 
results  reported  in  this  work  were  for  an  applied  load  which  was  oriented  normal  to 
the  original  crack,  mainly  for  simplicity  in  their  presentation.  That  is  not  a  limitation 
of  the  method,  except  that  the  crack  path  direction  predictor  is  still  controversial. 

There  are  several  interesting  results  that  come  from  the  parametric  studies 
shown  in  figures  in  the  previous  section.  First,  it  seems  that,  given  the  limited  number 
of  cases  examined  in  the  figures,  longer,  flatter  inclusions  may  tend  to  attract  cracks, 
whereas  shorter,  rounder  inclusions  may  tend  to  repel  cracks.  This  result  may  be 
significant  in  the  design  of  reinforcements  in  materials  if  the  energy  absorption  due  to 
crack  blunting  and  crack  path  deflection  are  intended  to  be  toughening  mechanisms. 


81 


Real  materials,  however,  have  more  than  one  inclusion,  and  a  naturally  growing  crack 
will  interact  with  all  of  the  inclusions  present.  The  method  presented  in  this  work, 
however,  can  be  extended  to  deal  with  multiple  inclusions  and/or  multiple  cracks.  The 
method  is  not  intended  to  provide  stochastic  information  about  cracks  in  reinforced 
materials.  Rather,  it  is  a  fundamental  building  block  of  what  could  become  a 
competent  micromechanics  model  for  dealing  with  cracks  near  inclusions  in  reinforced 
materials,  or  precipitation  hardened  metals. 

6.2  Limitations  of  Method 

The  inclusion  is  at  present  rigid  and  rigidly  bonded  to  the  surrounding 
medium.  In  the  predictions  given,  the  computation  was  forced  to  stop  if  the  crack 
attempted  to  enter  the  inclusion.  Also,  no  attention  was  paid  to  the  bond  strength  of 
the  inclusion  to  the  surrounding  medium.  In  real  materials,  when  a  crack  propagates 
to  an  interface  with  a  reinforcing  panicle,  one  of  three  things  can  happen.  Either  the 
panicle  can  debond  with  the  material,  or  the  crack  will  be  arrested  at  the  panicle,  or 
the  particle  itself  will  fracture.  Real  composite  materials  and  strengthened  metals 
exhibit  all  three  of  these  phenomena  during  fracture  processes.  The  methodology 
proposed  in  this  work  does  not  take  these  effects  into  account,  and  must  be  extended 
to  include  those  effects  if  that  is  deemed  to  be  critical  for  a  certain  application.  The 
results  shown  for  energy  absorption  for  cracks  near  open  holes  compared  to  cracks 
near  inclusions  is  evidence  of  this  limitation.  It  has  been  known  for  some  time  that 


drilling  a  hole  slightly  ahead  of  an  advancing  crack  tip  will  often  arrest  the  progression 
of  the  crack.  In  that  case  it  would  appear  that  the  hole  had  absorbed  energy  when  in 
fact  what  has  been  done  is  to  blunt  the  sharp  tip  of  the  crack,  or  in  the  terms  of  the 
work  presented  here,  to  remove  the  singularity  at  the  advancing  crack  tip. 

There  is  also  at  present  only  one  inclusion  included  in  the  analysis.  Most 
real  materials  of  interest  have  a  multitude  of  inclusions  or  reinforcing  particles  and/or 
fibers  in  a  multitude  ^  orientations.  The  method  can  be  extended  to  include  the  effect 
of  more  than  one  inclusion,  but  a  large  number  of  inclusions  would  be  difficult  to  treat 
with  this  method.  One  would  then  have  to  resort  to  stochastic  models  of  inclusion  or 
reinforcement  orientation  and  density  to  determine  the  energy  absorbed  by  crack  path 
deflection.  These  models  are  inherently  much  more  approximate  in  nature,  and  have 
yet  to  yield  real  insight  into  the  process  of  fracture  in  very  complicated  stress  fields. 

The  method  is  also  limited  to  regions  whose  stresses  are  two  dimensional 
in  nature  or  that  can  be  modelled  as  two  dimensional,  and  to  cases  in  which  the  two 
crack  faces  do  not  touch  or  try  to  overlap.  Potential  methods  are  usually  limited  to 
two  dimensions  as  the  potentials  are  significantly  easier  to  calculate.  Three 
dimensional  potential  methods,  while  they  exist,  are  not  commonly  used  because  of 
their  complexity.  If  the  stress  field  is  truly  three  dimensional,  the  techniques  described 
in  this  work  should  not  be  used.  Also,  if  the  crack  is  under  compression  in  that  its 
faces  would  touch,  an  extension  of  this  work  would  be  required. 
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6.3  Future  Areas  of  Investigation 

The  two  most  likely  areas  that  should  be  investigated  fc. lowing  this  work 
are  those  identified  as  limitations  to  the  method.  First,  for  those  inclusions  which  are 
relatively  flat,  with  cracks  that  are  attracted  to  them,  modelling  the  crack  propagation 
will  require  following  the  path  of  the  crack  either  through  the  inclusion  or  along  the 
boundary  between  the  inclusion  and  the  surrounding  material.  Each  of  these  cases 
would  have  to  be  treated  separately,  and  new  solutions  would  have  to  be  found  for  the 
stress  analysis  of  the  crack,  using  Green’s  Functions,  or  some  other  appropriate  metnod 
to  calculate  interaction  stresses.  Once  these  solutions  were  found,  they  could  be 
relatively  easily  implemented  within  the  structure  of  the  prediction  model  proposed  in 
this  work. 

The  second  area  of  future  investigation  that  seems  to  come  directly  from 
this  analysis  is  that  real  materials  that  are  reinforced  do  not  have  just  one  inclusion. 
The  potential  for  adding  multiple  inclusions  in  the  analysis  exists  within  the  framework 
of  the  analysis  presented  here.  Each  inclusion  would  merely  add  its  own  set  of  non¬ 
singular  terms  to  the  kernel  of  the  analysis,  and  each  would  also  add  stresses  on  the 
right  hand  sides  of  the  integral  equations.  Two  or  three  inclusions  would  be  a  fairly 
simple  task  to  handle  within  the  scope  of  the  technique  as  it  exists  now.  A  hundred 
or  more  inclusions  would  require  the  use  of  stochastic  techniques  which  are  outside 
the  scope  of  what  has  been  presented.  These  stochastic  techniques  could,  however,  be 
built  from  knowledge  gained  from  parametric  studies  of  a  crack  interacting  with  many 


fewer  inclusions.  In  any  case,  these  stochastic  models  would  have  to  be  much  more 
approximate  in  nature,  and  are  outside  the  present  scope.  Another  area  which  should 
be  investigated  is  some  type  of  self  consistent  micromechanics  model  of  multiple 
inclusions  and  multiple  cracks  in  an  elastic  material.  Some  of  these  techniques  are 
very  powerful,  and  should  be  explored  in  relation  to  the  solution  presented  here. 
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APPENDIX  A 


CRACK  PATH  PREDICTION  PROGRAM 


The  following  pages  contain  the  computer  program  used  to  perform  the 
crack  path  predictions  presented  in  this  work.  The  program  is  written  in  FORTRAN, 
and  can  be  compiled  as  is  to  function  on  most  computers. 
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on  noon  onn 


PROGRAM  CRACKY 

£  **************************************************************** 

C  THIS  PROGRAM  PREDICTS  THE  PATH  OF  A  NATURALLY  GROWING  CRACK 
C  WHICH  IS  INTERACTING  WITH  AN  ELLIPTIC  INCLUSION.  THE 
C  INCLUSION  IS  RIGID  AND  RIGIDLY  BONDED. 

Q  **************************************************************** 

c 

c 

Q  *************************************** ****•*••••***• *********** 

c 

C  WRITTEN  BY  EDWARD  M.  PATTON,  1987-1991 
C  IN  PARTIAL  FULFILLMENT  OF  THE  REQUIREMENTS  FOR 
C  THE  DOCTOR  OF  PHILOSOPHY  DEGREE  IN  MECHANICAL 
C  ENGINEERING  AT  THE  UNIVERSITY  OF  DELAWARE 
C 

C  HOPEFULLY  THE  VARIABLE  NAMES  THAT  I  HAVE  CHOSEN  WILL 
C  BE  APPARENT  TO  THE  READER  OF  THE  CODE.  IF  NOT,  THEN 
C  CAREFUL  READING  OF  MY  THESIS  SHOULD  ELUCIDATE. 

C 

Q  **************************************************************** 

c 

REALM  M.KAP 

COMPLEX  RKHAT(70),SIHAT(70),Zl,Z2,FACl,FAC2,zstr 
COMPLEX  XS  HE AR(70,70),YSHEAR (70,70) 

COMMON  /SPLIN/  X(70),Y(70),Y1(70),Y2(70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1,Z2,STRESS  J>HIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /RUSS/  YSTAR(70,70),YPSTAR(70,70),CAPY(70,70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(I40,]40) 
COMMON  /BLOK3/  SI(70),RK(70),N2,H,H2,N2P2,FAC1  ,FAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ,FJ,GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 

OPEN  NECESSARY  FILES  FOR  OUTPUT 

OPEN  (UNIT=l,FILE=’crackout’,FORM=’FORMATTED’) 

OPEN  (UNIT=2,file=’crackpaih’,FORM=’FORMATTED’) 

REWIND  1 
REWIND  2 

WRITE  HEADER  FOR  TECPLOT  -  YOU  CAN  USE  ANY  PLOTTING  PROGRAM 
THAT  YOU  LIKE  FOR  OFF-LINE  PLOTTING  OF  THE  CRACK  PATH. 

WRITE  (2,*)  ’TITLE  = . 

WRITE  (2,*)  ’VARIABLES  =  X  Y  ’ 

WRITE  (2,*)  ’ZONE’ 

SET  CONSTANTS 
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c 

R  =  1. 

PI  =  3.1415926536 
I  TIP  =  0 
ITRAP  =  0 
THETA  =  0. 

C 

C  BEGIN  THE  PROGRAM. 

C 

C  THE  FOLLOWING  SEQUENCE  OF  SUBROUTINE  CALLS 
C  SETS  UP  AND  SOLVES  THE  BOUNDARY  VALUE  PROBLEM  OF  A  CRACK 
C  INTERACTING  WITH  A  RIGID  ELLIPTICAL  INCLUSION. 

C 

C  FIRST  WE  GET  THE  REQUIRED  INPUT,  LIKE  WHERE  THE  CRACK  IS,  AND 
C  WHAT  SHAPE  THE  ELLIPSE  IS 
C 

CALL  INPUT(N) 

C 

C  NOW  THE  MAIN  LOOP  OF  THE  PROGRAM 

C  FIRST  WE  NEED  TO  PARAMETERIZE  THE  CRACK  ITSELF 

C 

100  CALL  EXTEND(N, ITRAP) 

C 

C  NOW  TO  SET  UP  THE  CONSTANTS  FOR  GERASOULIS’  METHOD 
C 

CALL  POINTS 
C 

C  BUILD  THE  MATRIX  AND  CALL  THE  MATRIX  INVERTER 
C 

CALL  BUILD(NB,IER) 

C 

C  IF  BUILD  OR  THE  MATRIX  INVERTER  SEES  AN  ERROR,  WE  STOP 
C 

IF  (IER.EQ.1)  GO  TO  9999 
C 

C  NOW  TO  CALCULATE  THE  K1  (PI  AND  P2)  AND  K2  (BT  AND  BT2) 

C  FOR  THE  CRACK.  THE  1  END  IS  THE  LEFT  END,  AND  THE 

C  2  END  IS  THE  RIGHT  END.  THE  CRACK  IS  ONLY  ALLOWED  TO 

C  GROW  FROM  THE  LEFT  END.  IF  YOPU  WANT  TO  GROW  FROM  THE 
C  RIGHT  END,  YOU  MUST  EDIT  SUBROUTINE  EXTEND 
C 

PI  =-COS(PHI(2))*RSS(l)*2.*PI/ABSFAC 
!  +SIN(PHI(2))*RSS(N2P2)*2.*P1/ABSFAC 

BT  =  (COS(PHI(2))*RSS(N2P2)+SIN(PHI(2))*RSS(l)) 

!  *2  *PI/ABSFAC 

P2  =-COS(PHI(N2P2))*RSS(n2p2- 1  )*2.*PI/ABSFAC 
!  +SIN(PHI(N2P2))*RSS(2*N2P2-2)*2.*PI/ABSFAC 

BT2  =  (COS(PHI(N2P2))*RSS(2*N2P2-2)+SIN(PHI(N2P2))*RSS(N2P2-l)) 

!  *2.*PI/ABSFAC 
C 
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NOW  CALCULATE  THE  ENERGY,  EVEN  THOUGH  WE  DO  NOT  USE  IT 

G  =  P1**2+BT**2 

WRITE  THE  RESULTS 

WRITE  (*,*)  Z 1  ,P  1 3T,P2,BT2 ,G ,AB S F AC ,N 2 
WRITE  (1  *)  Z1,P1,BT,P2,BT2,ABSFAC>J2 

CALL  THE  OUTPUT  ROUTINE 

CALL  OUT 

WRITE  COORDINATES  OF  NEW  CRACK  TIP  ON  RLE  2 

XTIP  =  REAL(SIHAT(2)) 

YTIP  =  AIMAG(SIHAT(2)) 

WRITE  (2,10)  XTIP,YTIP 
10  FORMAT  (G  15.6, 3X,G  15.6) 

NOW  FOR  THE  CRACK  PATH  PREDICTION.  THIS  IS  AN  ITERATIVE 
PROCESS  THAT  KEEPS  GOING  UNTIL  N  GETS  TOO  LARGE,  OR  THE 
CRACK  GROWS  PAST  THE  POINT  X  =  -3 

GROW  OFF  THE  LEFT  TIP 


ITIP  =  1 


TEST  TO  SEE  WHETHER  WE  HAVE  REACHED  OUR  LIMIT 
IF  (N2.LE.150.AND.REAL(Zl).GT.-3.)  THEN 
SET  CONSTANTS  AND  COUNTERS  FOR  THE  PREDICTION 


THOLD  =  0. 

BTOLD  =  .2 
IBACK  =  0 
THETA  =  0. 

DTHETA  =  4.*PI/180. 
NITTER  =  0 


CALCULATE  THE  POSITION  OF  THE  NEW  CRACK  TIP,  AND  SEE  IF  IT 
IS  INSIDE  THE  ELLIPSE 

110  ZSTR  =  Zl-CMPLX(COS(THETA)*DZ,SIN(THETA)*DZ) 

XTIP  =  REAL(ZSTR) 

YTIP  =  AIMAG(ZSTR) 

ZSANG  =  ATAN(YTIP/XTIP) 

ZSTRMAG  =  SQRT(XTIP*  *2+  YTIP*  *2) 

RADIUS  =  SQRT((L-M*M)**2/(COS(ZSANG)**2*(L-M)**2 
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'  _.IN(ZSANG)**2*(1.+M)**2)) 

IF  THE  CRACK  TIP  IS  IN  THE  ELLIPSE,  STOP 

IF  (ZSTRMAG-RADIUS.LT.O.)  GO  TO  9999 

OTHERWISE,  GO  ON  TO  THE  ITERATIVE  SOLUTION.  I  HAVE 
GIVEN  THIS  THING  AN  ARBITRARY  MAXIMUM  OF  50  ITERATIONS 
TO  FIND  WHERE  IT  WANTS  TO  GO.  IF  THE  PROBLEM  IS  POSED 
CORRECTLY,  YOU  SHOULD  NEED  NO  MORE  THAN  15  ITERATIONS, 
WITH  AN  AVERAGE  OF  ABOUT  10  PER  INCREMENT  OF  GROWTH 

DO  200  ITRAP  =  1,50 

SAME  SET  OF  SUBROUTINES.  WHAT  THIS  DOES  IS  TO  TRY 
SEVERAL  DIFFERENT  DIRECTIONS  UNTIL  IT  FINDS  THE 
DIRECTION  OF  MINIMUM  SHEAR  STRESS.  THEN  IT  GOES 
BACK  TO  STATEMENT  100  AND  STARTS  FRESH  WITH  A  NEW 
CRACK  GROWTH  INCREMENT 

CALL  EXTEND(N, ITRAP) 

CALL  POINTS 

CALL  BUILD(NB.IER) 

IF  (IER.EQ.l)  GO  TO  9999 

PI  =-COS(PHI(2))*  RSS(  I  )*2.*PI/ABSFAC 
!  +SIN(PHI(2))*RSS(N2P2)*2.*PI/ABSFAC 

bt  =  (COS(PHI(2))*RSS(N2P2)+SIN(PHl(2))*RSS(l)) 

!  *2.*PI/ABSFAC 
G  =  P1**2+BT**2 

LET  THE  USER  KNOW  THAT  YOU  HAVE  COMPLETED  AN  INTERATION, 
AND  SHOW  HIM  THE  PERTINENT  RESULTS 

WRITE  (*,*)  P  1,BT,G, THETA*  180./PI 

NOW,  IF  THE  NEW  K2  IS  LESS  THAN  THE  OLD  K2,  DIVIDE  DTHETA 
IN  HALF  AND  GO  BACKWARDS.  IF  YOU  HAVEN’T  BEEN  BACKWARDS 
TWICE  YET,  YOU  HAVEN’T  FOUND  THETA 

IF  (ABS(BT).GT.BTOLD.AND.IBACK.LT.2)  THEN 
DTHETA  =  -DTHETA/2. 

THOLD  =  THETA 

j  THETA  =  THETA+DTHETA 

PUT  SOME  REASONABLE  LIMITS  ON  THETA 
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IF  (THETA*  180./PI.GT.90.)  THEN 
THETA  =  90  *PI/180. 

GO  TO  201 

ELSEIF(THETA*  180./PI.LT.-90.)  THEN 
THETA  =  -90.*PI/180. 

GO  TO  201 
ENDIF 

CHECK  TO  SEE  IF  YOU  ARE  INSIDE  THE  ELLIPSE 

ZSTR  =  SIHAT(4)-CMPLX(COS(THETA)*DZ,SIN(THETA)*DZ) 
XTIP  =  REAL(ZSTR) 

YTIP  =  AIMAG(ZSTR) 

ZSANG  =  ATAN(YTIP/XTIP) 

ZSTR  MAG  =  SQRT(XTIP**2+YTIP**2) 

RADIUS  =  SQRT((l.-M*M)**2/(COS(ZSANG)**2*(l.-M)**2 
+SIN(ZSANG)**2*(1.+M)**2)) 

IF  (ZSTRMAG-RADIUS.LT.C.)  GO  TO  9999 

YOU  JUST  WENT  BACKWARDS,  SO  INCREMENT  IBACK 

IBACK  =  IBACK+1 

RESET  THE  VALUE  OF  OLD  K2 

BTOLD  =  ABS(BT) 

AND  COUNT  THE  ITERATION 

NITTER  =  NITTER+1 

TRY  AGAIN 

GO  TO  200 

OR  IF  YOU  ARE  STILL  GOING  DOWN  IN  MAGNITUDE,  GO  THE  SAME 
DIRECTION  AND  RESET  OLD  K2 

ELSEIF  (ABS(BT).LE.BTOLD)  THEN 
BTOLD  =  ABS(BT) 

THOLD  =  THETA 

THETA  =  THETA+DTHETA 

SET  LIMITS  ON  THETA 

IF  (THETA*  180./PI.GT.90.)  THEN 
THETA  =  90.*  PI/ 180. 

GO  TO  201 

ELSEIF(THETA*180./PI.LT.-90.)  THEN 
THETA  =  -90.*PI/180. 
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GO  TO  201 
ENDIF 
C 

C  CHECK  TO  SEE  IF  YOU  HAVE  GONE  IN  THE  ELLIPSE 
C 

ZSTR  =  SIHAT(4)-CMPLX(COS(THETA)*DZ,SIN(THETA)*DZ) 

XT  IP  =  REAL(ZSTR) 

YTIP  =  AIMAG(ZSTR) 

ZSANG  =  ATAN(YTIP/XTIP) 

ZSTR  MAG  =  SQRT(XTIP**2+YTTP**2) 

RADIUS  =  SQRT((l.-M*M)**2/(COS(ZSANG)**2*(l.-M)**2 
!  +SIN(ZSANG)**2*(1.+M)**2)) 

IF  (ABS(RADIUS-ZSTRMAG).LT. .00001 )  GO  TO  9999 
C 

C  COUNT  THE  ITERATION  AND  GO  AGAIN 
C 

NITTER  =  NITTER+1 
GO  TO  200 
C 

C  OR,  IF  THE  NEW  K2  IS  BIGGER,  AND  YOU  HAVE  BEEN  BACKWARDS 
C  TWICE,  YOU  HAVE  FOUND  THETA.  INCREMENT  THE  CRACK  AND 
C  GO  OUT  OF  THE  INNER  LOOP 
C 

ELSEIF  (ABS(BT).GT.BTOLD.AND.IBACK.EQ.2)  THEN 
THETA  =  THOLD 
GO  TO  201 
ENDIF 
C 

CALL  OUT 
C 

200  CONTINUE 
C 

201  CONTINUE 
C 

C  TELL  THE  USER 
C 

WRITE  (*,*)  ’AT  N*2  =  ’.N2,’  THETA  =  ’ .THETA*  180./PI, 

!  ’  ITERATIONS  ’.NITTER 
C 

C  AND  GO  BACK  TO  THE  BEGINNING 
C 

GO  TO  100 
C 

C  HERE  IS  WHERE  YOU  GO  IF  N  IS  TOO  BIG  OR  IF  THE  CRACK 
C  IS  FAR  ENOUGH  TO  THE  LEFT 
C 

ENDIF 

C 

C  ALL  DONE 
C 
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9999  CONTINUE 
CLOSE  (UNIT=I) 
CLOSE  (UNIT=2) 
STOP 
END 


SUBROUTINE  INPUT(N) 

***************************************************************** 

THIS  ROUTINE  PROVIDES  ALL  OF  THE  REQUIRED  INPUT  TO  CRACKY 
********************************* ************** ******* *********** 


REALM  MJKAP 

COMPLEX  RKHAT(70),SIHAT(70),Z  1  JZ2.FAC  1  ,FAC2 
COMPLEX  XS  HEAR(70,70), YSHE AR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1,Z2,STRESS,PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  Sl(70),RK(70)(N2,H,H2,N2P2,FACl,FAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ.FJ.GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 

CONVER  =  3.1415926536/180. 

200  FORMAT  (2G10.4) 

CRACK  HAS  A  LENGTH  OF  1  UNIT,  WITH  LEFT  TIP  AT  Z1  AND 
RIGHT  TIP  AT  Z2 


CLEN  =  1. 

WRITE  (*,*)  ’  WHERE  IS  THE  CLOSE  TIP  OF  THE  CRACK  ’ 

READ  (*,200)  Z1 

WRITE  (*,*)  ’  NOW  WHAT  IS  THE  ELLIPTICITY  RATIO  OF  THE  ELIPSE  (0  TO  1)’ 
READ  (*,200)  M 

WRITE  (*,*)  ’  WHAT  IS  THE  ANGLE  OF  THE  CRACK  ’ 

READ  (*,200)  PHIIN 
PHI  IN  =  PHIIN*CONVER 
SINPHI  =  SIN(PHIIN) 

COSPHI  =  C0S(PH1IN) 

SET  UP  THE  CRACK 


Z2  =  C  MPLX  (RE  AL(Z  1  )+CLEN*COS  PHI  ,AI  MAG(Z  1  )+CLEN*  SINPHI) 

WRITE  (*,100)  Z1 .Z2.PHIIN/CONVER 
100  FORMAT  (’  Z1  =  ’.2G10.4,’  Z2  =  ’.2G10.4,’  PHI  =  \G10.4) 

C 

C  UNCOMMENT  THE  COMMENTED  LINES  AND  COMMENT  OUT  THE  ASSIGNMENTS 
C  IF  YOU  WANT  TO  HAVE  MORE  CONTROL  OVER  THE  INPUT 
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c 

C  WRITE  (*  ,*)  ’  WHAT  IS  THE  ANGLE  BETWEEN  THE  CRACK  AND  THE  LOAD’ 

C  READ  (*,200)  ALPHA 
C  ALPHA  =  ALPHA*CONVER+PHIIN 
ALPHA  =  90.*CONVER+PHIIN 
N  =  5 

C  WRITE  (* ,*)  ’  WHAT  WOULD  YOU  LIKE  FOR  THE  VALUE  OF  N’ 

C  READ  (*,*)  N 

WRITE  (* ,*)  ’  WHAT  WOULD  YOU  LIKE  FOR  KAPPA’ 

READ  (*,200)  KAP 
C  KAP  =  1.67 
C 

C  NOW  SET  UP  THE  CRACK  PARAMETERS  THAT  WE  NEED  FOR  THE  SOLUTION 
C 

FAC1  =  (Z2+Zl)/2. 

FAC2  =  (Z2-Zl)/2. 

ABSFAC  =  SQRT  (RE  AL(FAC2)*  *  2+ AIMAG(F  AC2)*  *2) 

DZ  =  ABSFAC/FLOAT(N) 

C 

C  DONE 
C 

RETURN 

END 

C 

SUBROUTINE  OUT 
C 

c 

REAL*4  M,KAP 

COMPLEX  RKHAT(70),SIHAT(70),Z1  ,Z2,FAC  1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1 ,Z2, STRESS  .PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZJ>HI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 

COMMON  /BLOK3/  SI(70),RK(70)^2,HJ32,N2P2,FAC1  JFAC2, ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ,FJ,GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 
C 

Q  *************************************************************** 

C  THIS  ROUTINE  PERFORMS  ALL  OF  THE  OUTPUT  NECESSARY  FOR  CRACKY 
£  *************************************************************** 

C 

C 

WRITE  (1,*)  ’  Output  from  Crack  Program’ 

C 

WRITE  (1,*) 

C 

WRITE  (1,200)  (RSS(I)*2.*3.1415926536/ABSFAC,I=l,2*N2P2-2) 

WRITE  (1,*) 

99 


nnnonn  oonno  nnoonononon  non 


WRITE  (1,200)  (RHS (I),I=  1 ,2*N2P2-2) 
200  FORMAT  (G15.5) 

C 

WRITE  (1  *) 

DONE 

RETURN 
END 


SUBROUTINE  EXTEND (N.ITRAP) 


********************************************************** 

THIS  ROUTINE  PROVIDES  THE  PROPER  VALUES  OF  THE  COORDINATES 
OF  THE  EXTENDED  CRACK  IN  TABULAR  FORM,  AND  ALSO  INCREASES 
N  BY  THE  APPROPRIATE  AMOUNT  (1),  AT  AN  ANGLE  CALCULATED 
IN  FUNCTION  CRACK 

********************************************************** 


REAL*4  M.KAP 

COMPLEX  RKHAT(70),SIHAT(70),Z1,Z2,FAC1,FAC2 
COMPLEX  XS  HE AR(70,70), Y S  HE AR (70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1,Z2,STRESS,PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /SPLIN/  X(70),Y(70),Y1(70),Y2(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70)^2,H,H2,N2P2,FAC1  JFAC2.ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ,FJ,GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 

ONLY  INCREMENT  N  IF  WE  HAVE  BEEN  THROUGH  THE  FIRST  SOLUTION 
AND  AT  THE  FIRST  ITERATION  LOOKING  FOR  THETA  (ITRAP  =  1) 


IF  (ITIP.EQ.  1  .AND.ITRAP.EQ.  1)  N=N+1 

N2  =  N*2 
N2P2  =  N2+2 
N2P2M1  =  N2P2-1 
H  =  1VFL0AT(N) 

FOR  ITERATIONS  AFTER  THE  FIRST,  BACK  OFF  TO  WHERE  THE  CRACK 
TIP  WAS  AT  THE  START  OF  THIS  INCREMENT 

IF(ITRAP.GT.l)  Z1  =  SIHAT(4) 

FIRST  TIME  THROUGH  FOR  THE  WHOLE  PROGRAM,  SET  UP  THE  CRACK 
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c 


IF  (ITIP.EQ.0)  THEN 


DO  1  I  =  2.N2P2 

SIHAT(I)  =  ((FL0AT(I)-2.)*H-1.)*FAC2+FAC1 

IF  (I.EQ.2)  SIHAT(I)  =  Z1 

IF  (I.EQ.N2P2)  SIHAT(I)  =  Z2 

PHI(I)  =  PHI  IN 

X(I-1)  =  REAL(SIHAT(I)) 

Y(I-l)  =  AIMAG(SIHAT(I)) 

Yia-1)  =  TAN(PHI(I)) 

1  CONTINUE 

CALL  THE  SPLINE  ROUTINE 

CALL  SPLINE  (N2P2M1,X,Y,Y1,Y2) 

DO  2  K  =  2.N2P2 

RKHAT(K)  =  (((FL0AT(K)-2.)*H-1.)+H/2.)*FAC2+FAC1 

2  CONTINUE 

FOR  SUBSEQUENT  TRIPS  THROUGH  THE  PROGRAM,  WE  NEED  TO 
ADD  AN  INCREMENT  IN  THE  SPECIFIED  DIRECTION,  AND  THEN 
GET  THE  FIRST  AND  SECOND  SPATIAL  DERIVATIVES  OF  THE 
CRACK  BY  CALLING  SPLINE 

ELSEIF  (ITIP.EQ.l)  THEN 

FIRST  ITERATION  AGAIN,  SHOVE  ALL  OF  THE  CRACK  GEOMETRY 
TWO  SPACES  HIGHER  IN  EACH  OF  THE  ARRAYS  THAT  HOLD  CRACK 
GEOMETRY 

IF  (ITRAP.EQ.l)  THEN 
DO  5  J=2JM2 
I  =  N2P2+2-J 
RKHAT(I)  =  RKHAT(I-2) 

SIHATfl)  =  SIHAT(I-2) 

PHI(I)  =  PHia-2) 

X(I-1)  =  X(I-3) 

Y(I-l)  =  Y(I-3) 

Yl(I-l)  =  Y  1(1-3) 

5  CONTINUE 
ABSFAC  =  ABSFAC+DZ 
ENDIF 

THIS  PART  ACTUALLY  PERFORMS  THE  INCREMENT 

PHI(3)  =  -THETA 
PHI(2)  =  -THETA 

SIHAT(2)  =  Zl-CMPLX(COS(PHI(2))*2.*DZ,SIN(PHI(2))*2.*DZ) 

SIHAT(3)  =  Zl-CMPLX(COS(PHI(3))*DZ,SIN(PHI(3))*DZ) 
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RKHAT(2)  =  Zl-CMPLX(COS(PHI(2))*3.*D272.,SIN(PHI(2))*3.*DZ/2.) 
RKHAT(3)  =  Zl-CMPLX(COS(PHI(3))*DZ/2.,SIN(PHI(3))*DZ/2.) 

X(l)  =  REAL(SIHAT(2)) 

X(2)  =  REAL(SIHAT(3)) 

Y(l)  =  AIMAG(SIHAT(2)) 

Y(2)  =  AIMAG(SIHAT(3)) 

Yl(l)  =  TAN  (PHI  (2)) 

Yl(2)  =  TAN(PHI(3)) 

Z1  =  SIHAT(2) 

GET  OUR  SPATIAL  DERIVATIVES 

CALL  SPLINE  (N2P2M1,X,Y,Y1,Y2) 

C 

ENDIF 

C 

RETURN 

END 

C 

SUBROUTINE  SPLINE  (N,X,Y,Y1,Y2) 

C 

Q  ************************************************************* 

C  THIS  ROUTINE,  EXERPTED  FROM  "NUMERICAL  RECIPES  IN  FORTRAN", 
C  CALCULATES  THE  LOCAL  FIRST  AND  SECOND  DERIVATIVES  FOR  USE 
C  IN  SUBROUTINE  RUSSIAN 

( 2  ************************************************************* 

C 

DIMENSION  U(70),X(70),Y (70), Y 1  (70),Y2(70) 

C 

YP1  =  l.e35 
YPN  = 1x35 
DO  10  IN  =  1,N 
Y2(IN)  =  0. 

10  CONTINUE 

IF  (YP1.GT..99E30)  THEN 
Y2(l)  =  0. 

U(l)  =  0. 

ELSE 

Y2(l)  =  -.5 

U(l)  =  (3  ,/(X(2)-X(  1  )))*((  Y(2)- Y(  1  ))/(X(2)-X(  1 ))- YP 1 ) 

ENDIF 

DO  11  I=2,N-1 

SIG  *  (X(I)-X(I- 1  ))/(X(I+ 1)-X(I- 1)) 

P  =  SIG*Y2(I-l)+2. 

Y2(I)  =  (SIG-1.)/P 

U(I)  =  (6.*((Y  (1+ 1)- Y  (I))/(X(I+ 1  )-X(I))-(  Y  (I)- Y  (I- 1 )) 

!  /(X(I)-X(I- 1  )))/(X(I+ 1  )-X(I-l  ))-SIG*U(I-  1))/P 

11  CONTINUE 
IF(YPN.GT.  .99E30)  THEN 

QN  =  0. 
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UN  =  0. 

ELSE 
QN  =  .5 

UN  =  (3./(X(N)-X(N-l)))*(YPN-(Y(N)-Y(N-l))/(X(N)-X(N-l))) 
END  IF 

Y2(N)  =  (UN-QN*U(N-1))/(QN*Y2(N-1)+1.) 

IX)  12  K  =  N-1,1,-1 

Y2(K)  =  Y2(K)*Y2(K+1)+U(K) 

12  CONTINUE 
RETURN 
END 


SUBROUTINE  POINTS 

********************************************* ********* ********* 
CALCULATES  THE  SI  AND  RK  QUANTITIES  FOR  GERASOULIS’ 
TECHNIQUE.  ALSO  SETS  UP  THE  WEIGHT  FUNCTIONS  BY  CALLING 
SUBROUTINE  VWGT  AND  WTFUNS 
*************************************************************** 


REALM  M.KAP 

COMPLEX  RKHAT(70),SIHAT(70),Z1,Z2,FAC  1 TAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1(Z2,STRESS,PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70)JM2,H,H2,N2P2,FACU;AC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ,FJ,GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 

N2P2MI  =  N2P2-1 

CALCULATE  SI  ARRAY 

DO  1  I  =  2.N2P2 

SI(I)  =  (FLOAT(I)-2.)*H- 1 . 

SI(2)  =  -1. 

SI(N2P2)  =  1. 

1  CONTINUE 

CALCULATE  RK  ARRAY 

DO  2  K  =  2N2P2 
RK(K)  =  SI(K)+H/2 

2  CONTINUE 

SET  CONSTANT  WEIGHTS 
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c 


CALL  VWGT 


C 

DO  3  K  =  2JM2P2M1 
DO  4  I  =  2.N2P2 

CALL  WTFUNS(RK(K)J) 
WnC(K,I)  =  EJ+FJ+GJ 
4  CONTINUE 
3  CONTINUE 
C 

RETURN 

END 


SUBROUTINE  WTFUNS(S.K) 

ft************************************************************** 

THIS  SUBROUTINE  COMPUTES  THE  APPROPRIATE  WEIGHT  FUNCTIONS  FOR 
THE  PIECE-WISE  SMOOTH  POLYNOMIAL  APPROXIMATION 
*************************************************************** 

REAL*4  M,KAP 

COMPLEX  RKHAT(70),SIHAT(70),Z  1.Z2.FAC  1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHAyZl  ,Z2, STRESS  .PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 

COMMON  /BLOK3/  SI(70),RK(70),N2,H,H2,N2P2,FAC1,FAC2)ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ,FJ,GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 


CHECK  TO  SEE  IF  K  IS  ODD  OR  EVEN 

KCHEK1  =  K/2 
KCHEK2  =  2*KCHEK1 
DF  (KCHEK2.NE.K)  THEN 

K  IS  ODD 

0  =  0. 

GJ  =  0. 

J  =  (K+l)/2 

FJ  =  -RM(S,SI(2*J),SI(2*J-2)J)/H2 
RETURN 
C 

ELSEIF  (K.EQ.2)  THEN 
C 
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K  =  2 

FJ  =  0. 

GJ  =  0. 

J  =  (K+2)/2 

EJ  =  RM(S,SI(2*J),SI(2*J-1)J)/(2.*H2) 
RETURN 

ELSEIF  (K.EQ.N2P2)  THEN 

K  =  N2P2 

FJ  =  0. 

EJ  =  0. 

J  =  K/2 

GJ  =  RM(S,SI(J*2-1),SI(J*2-2)J)/(2.*H2) 
RETURN 

ELSE 

K  IS  EVEN 
J  =  (K+2)/2 

EJ  =  RM(S,SI(J*2),SI(J*2-1)J)/(2.*H2) 

J  =  K/2 

GJ  =  RM(S,SI(J*2-1),SI(J*2-2)J)/(2.*H2) 
FJ  =  0. 

RETURN 

ENDIF 

END 


SUBROUTINE  VWGT 
C 

Q  ************************************************************* 

C  THIS  ROUTINE  COMPUTES  THE  CONSTANTS  FOR  THE  POLYNOMIAL 
C  APPROXIMATION  -  SEE  GERASOULIS’  PAPER 

Q  ************************************************************* 

c 

REALM  M,KAP 

COMPLEX  RKHAT(70),SIHAT(70),Z 1  ,Z2pAC  1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHArZl  ^.STRESS  PHI  IN 
COMMON  /GEOM2/  mP,THETA,DZPHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70),N2,H,H2,N2P2,FAClpAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ.FJ.GJ 
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COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 


H2  =  H**2 
NP1  =  N2P2/2 

DO  1  J1  =  2N2P2 

PSI(Jl)  =  ASIN(SI(J1)) 

1  CONTINUE 

DO  2  J2=  1 JNP1 

TAU(J2)  =  TAN(PSI(J2*2)/2.) 

2  CONTINUE 

DO  3  J3  =  2JM2P2 
J3CHK1  =  J3/2 
J3CHK2  =  J3CHK1*2 

J3  IS  ODD 

IF  (J3.NE.J3CHK2)  THEN 

AJ  =  0. 

CJ  =  0. 

J3P  =  (J3+l)/2 

BJ  =  -HCON(J3P,SI(J3P*2),SI(J3P*2-2))/H2 
GO  TO  4 

J3  =  2 

ELSEIF  (J3.EQ.2)  THEN 
J3P  =  2 

AJ  =  HCON(J3P,SI(J3P*2),SI(J3P*2-l))/(2.*H2) 
BJ  =  0. 

CJ  =  0. 

GO  TO  4 

J3  =  N2P2 

ELSEIF  (J3.EQ.N2P2)  THEN 
J3P  =  NP1 

CJ  =  HCON(J3P,SI(J3P*2-l),SI(J3P*2-2))/(2.*H2) 
AJ  =  0. 

BJ  =  0. 

GO  TO  4 

J3  IS  EVEN 
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ELSE 

C 

BJ  =  0. 

J3P  =  (J3+2)/2 

AJ  =  HCON(J3P,SI(J3P*2),SI(J3P*2-l))/(2.*H2) 
J3P  =  J3/2 

O  =  HCON(J3P,SI(J3P*2-l),SI(J3P*2-2))/(2  *H2) 
C 

ENDIF 

C 

4  VI(J3)  =  AJ+BJ+a 
C 

3  CONTINUE 
RETURN 
END 


FUNCTION  HCON(J,X,Y) 

************************************************************* 
THIS  FUNCTION  IS  USED  TO  DETERMINE  THE  WEIGHTS 
FOR  THE  PIECEWISE  SMOOTH  POLYNOMIAL  APPROXIMATION 
************************************************************* 

REALM  M.KAP 

COMPLEX  R  KH  AT  (70)  ,S  I  HAT  (70)  JZ 1  ,Z2,FAC  1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Zl,Z2,STRESStPHlIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70)Nf2,H,H2,N2P2,FACl>FAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6 /  EJ.FJ.GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 


J2  =  J*2 

TERM1  =  (X*Y+.5)*(PSI(J2)-PSI(J2-2)) 

TERM2  =  (X+Y)*(COS(PSI(J2))-COS(PSI(J2-2))) 
TERM3  =  (SIN(2  *PSI(J2))-SIN(2.*PSI(J2-2)))/4. 
HCON  =  TERM  1+TERM2-TERM3 
RETURN 
END 


FUNCTION  RM(S,X,YJ) 

C 

Q  ************************************************************* 

C  THIS  FUNCTION  IS  CALLED  BY  POINTS  AND  IS  USED 


107 


on  n  n  n  on 


C  IN  DETERMINING  WEIGHTS  FOR  THE  POLYNOMIAL 
C  APPROXIMATION 

Q  ************************************************************* 

c 

REAL* 4  M,KAP 

COMPLEX  RKHAT(70),SIHAT(70),Z1  Z2.FAC1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1.Z2<STRESS  PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70)VN2,H,H2,N2P2,FACI^AC2,ABSFAC 
COMMON  /BLOK4/  PSI(70)  ,TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ.FJ.GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM{70,70),XSHEAR,YSHEAR 


JM1  =  J-l 
J2  =  J*2 
J2M2  =  J2-2 

AINT1  =  SQRT(1.-S**2) 

AINT2  =  (- 1 .+ AINT1 +S*T  AU(J))*TAU(JM  1  )+S-(  1 .+ AINT1  )*T  AU(J) 
AINT3  =  (-1.-AINT1+S*TAU(J))*TAU(JM1)+S-(1.-AINT1)*TAU(J) 
AINT4  s  ABS(AINT2/AINT3) 

AJOFS  =ALOG(AINT4)/AINTl 

TERM1  =  SQRT(1.-SI(J2)**2) 

TERM2  =  SQRT(1.-SI(J2M2)**2) 

TERM3  =  (S**2+X*Y-S*(X+Y))*  AJOFS 
TERM4  =  (S-X-Y)*(PSI(J2)-PSI(J2M2)) 

RM  =  -TERM1+TERM2+TERM3+TERM4 

RETURN 

END 


SUBROUTINE  BUILD  (NB.IER) 

C 

Q  ************************************************************* 

C  THIS  ROUTINE  BUILDS  THE  MATRIX  AND  PERFORMS  THE 
C  INVERSION  OF  THE  MATRIX  TO  SOLVE  THE  PROBLEM 
£  ************************************************************* 
C 

REALM  M.KAP 

COMPLEX  RKHAT(70),SIHAT(70)rZI22,FACl,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1,Z2,STRESSPHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /SPLIN/  X(70),Y(70),Y1(70),Y2(70) 
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COMMON  /RUSS/  YSTAR(70,70),YPSTAR(70,70),CAPY(70,70) 
COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140).ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70),N2,H,H2,N2P2,FAC1,FAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ.FJ.GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR>YSHEAR 
COMPLEX  SFORCEPFORCE 

ZERO  OUT  THE  MATRICES  AND  RIGHT  HAND  SIDES,  INCLUDING 
THE  RUSSIAN  STUFF 

DO  30  IA  =  1.N2P2 
RHS(IA)  =  0. 

RSS(IA)  =  0. 

DO  30  JA  =  1JM2P2 
ALHS(IAJA)  =  0. 

ALSS(IAJA)  =  0. 

YSTAR(IA3A)  =  0. 

YPSTAR(IAJA)  =  0. 

CAPY(IAJA)  =  0. 

CONTINUE 
N2P2M1  =  N2P2-1 
NB  =  2*N2P2M1 

NOW  SET  UP  THE  INTERACTION  TERMS  FROM  THE  KERNEL 


DO  50  I  =  2.N2P2 
DO  60  J  =  2J^2P2M1 
CALL  KERNEL(SIHAT(I)J?KHAT(J),I,J) 
CONTINUE 
CONTINUE 


AND  NOW  THE  STUFF  FROM  THE  RUSSIAN  PAPER 
CALL  RUSSIAN  (N2P2) 


NOW  LOOP  THROUGH  AND  SET  UP  THE  LEFT  AND  RIGHT  SIDES 
OF  THE  EQUATION  (A)X(SOLUTION)  =  (B) 

DO  1  JN  =  1JM2P2M1 
J  =  JN+1 
IT  =  JN 

GET  THE  LOCAL  ANGLES  RIGHT.  AND  THFJR  ctnes  AND  COSINES 

SINPHI  =  SIN(PH1(J)) 

COSPHI  =  COS(PHI(J)) 
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COS  3  PHI  =  COS(3.*PHI(J)) 

SIN3PHI  =  SIN(3.*PHI(J)) 

JT  =  JN+N2P2M1 

IX)  THE  RIGHT  HAND  SIDE  FIRST 

IF(J.EQ.N2P2)  THEN 
RHS(JN)  =  0. 

RHS(JT)  =  0. 

ELSE 

RHS(JN)  =  COSPHI**2*REAL(PFORCE(RKHAT(J)) 

!  +SFORCE(RKHAT(J)))+SINPHI**2*REAL(PFORCE(RKHAT(J)) 

!  -SFORCE(RKHAT(J)))-2.*SINPHI*COSPHI 
!  *AIMAG(SFORCE(RKHAT(J))) 

RHS(JT)  =  -2.*SlNPHI*COSPHI*(REAL(SFORCE(RKHAT(J)))) 

!  +(SINPHI**2-COSPHI**2)*AIMAG(SFORCE(RKHAT(J))) 

ENDIF 

NOW  THE  REST  OF  THE  MATRIX 

DO  1  IY  =  1JM2P2M1 
1  =  IY  +  1 
ITAU  =  IY 
IX  =  IY+N2P2MI 

FIRST  THE  KERNEL  STRESSES 

XSIGX  =  XNORM(IJ)-REAL(XSHEAR(U» 

XSIGY  =  XNORM(U)+REAL(XSHEAR(U)) 

XSIGXY  =  AIMAG(XSHEAR(U)) 

YSIGX  =  YNORM(U)-REAL(YSHEAR(U)) 

YSIGY  =  YNORM(IJ)+REAL(YSHEAR(U)) 

YSIGXY  =  AIMAG(YSHEAR(U)) 

NOW  THE  ANGULAR  COEFFICIENTS  FOR  AN  ARBITRARILY  ORIENTED 
STRAIGHT  CRACK 

BYX=  COSPHI+COS3PHI 
BYY  *  3*COSPHI-COS3PHI 
BYXY  =  SIN3PHI-SINPHI 
BXX  =  -(3*SINPHI+SIN3PHI) 

BXY  =  SIN3PHLSINPHI 
BXXY  =  COS3PHI+COSPHI 
SNY  =  SINPHI**2*BYX+COSPHI**2*BYY 
-2.*SINPHI*COSPHI*BYXY 
SNX  =  SINPHI**2*BXX+COSPHI**2*BXY 
!  -2*SINPHI*COSPHI*BXXY 

STY  =  SINPHI*COSPHI*(BYX-BYY) 

!  +BYXY*(SINPHI**2-COSPHI**2) 

STX  =  SINPHI*COSPHI*(BXX-BXY) 


no 
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!  +BXXY*(SINPHI**2-COSPHI**2) 

NOW  THE  RUSSIAN  STUFF  FOR  A  NON-STRAIGHT  CRACK 

BOTTOM  =  (L+YSTAR(IT,ITAU)**2)**2 

RNX  =  CAPY(IT,ITAU)*(YSTAR(IT,ITAU)**2-L-Y1(IT)* 

!  Y  STAR(IT  JTAU)*(3.+ Y  STAR(IT,rTAU)*  *2))/BOTTOM- 

!  YPSTAR(IT.ITAU) 

RNY  =  (YPSTAR(IT,ITAU)-CAPY(IT,ITAU))*(YSTAR(IT,ITAU)**2-1.) 

!  *(Y  STAR(IT  JTAU)- Y 1  (IT))/BOTTOM 

RTX  =  (YPSTAR(IT,ITAU)-CAPY(IT,ITAU))*(YSTAR(IT,ITAU) 

!  *(3.+YSTAR(IT,ITAU)**2)+Yl(IT)*(YSTAR(IT,ITAU)**2-l.)) 

!  /BOTTOM 

RTY  =  (YPSTAR(ITJTAU)-CAPY(IT,ITAU))*(YSTAR(IT,ITAU)**2-1.) 

!  *(  1 .+ Y 1  (IT)*  YSTAR(IT,ITAU))/BOTTOM+ YPSTAR(IT,IT  AU) 

NOW  BUILD  THE  MATRIX  ALHS  THAT  IS  TO  BE  SENT  OFF  TO  THE 
GAUSS -JORDAN  MARTI X  INVERSION  ROUTINE 

IF  (J.EQ.N2P2)  THEN 
ALHS(JN.IY)  =  VI(I)*COSPHI 
ALHS(JN.IX)  =  -VI(I)*SINPHI 
ALHS(JT,IY)  =  VI(I)*SINPHI 
ALHS(JT,IX)  =  VI(I)*COSPHI 
ELSE 

ALHS(JN.IY)  =  SNY*WIK(J,I)/ABSFAC-VI(I)*(YSIGX*SINPHI**2+ 

!  YSIGY*COSPHI**2-2.*YSIGXY*SINPHI*COSPHI) 

!  -VI(I)*RNY 

ALHS(JT,IY)  =  ST  Y*  WIK(J  ,1)/ ABSF  AC  -VI  (I)*  (SINPH1*  COSPH1* 

!  (YSIGX-YSIGY)+YSIGXY*(SINPHI**2-COSPHI**2)) 

!  -VI(I)*RTY 

ALHS(JN,IX)  =  SNX*WIK(J,I)/ABSFAC-VI(I) 

!  *(XSIGX*SINPHI**2+XSIGY*COSPHI**2 

!  -2.*XSIGXY*SINPHI*COSPHI) 

!  -VI(I)*RNX 

ALHS(JT,IX)  =  STX*WIK(J,I)/ABSFAC-VI(I) 

!  *(SINPHI*COSPHI*(XSIGX-XSIGY) 

!  +XSIGXY*(SINPHI**2-COSPHI**2)) 

!  -VI(I)*RTX 

ENDIF 

1  CONTINUE 

NOW  SAVE  THE  LEFT  HAND  SIDE  OF  THE  EQUATION,  AND 
SAVE  THE  FORCING  FUNCTIONS  FROM  THE  RIGHT  HAND  SIDE 

DO  110  I=1JMB 
RSS(I)  =  RHS(I) 

DO  110  K=1,NB 
ALSS(I.K)  =  ALHS(I,K) 

1 10  CONTINUE 
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c 

C  NOW  CALL  THE  MATRIX  INVERSION  ROUTINE  -  ALSS  GETS 
C  REPLACED  WITH  ITS  INVERSE,  AND  RSS  GETS  REPLACED 
C  WITH  THE  SOLUTION  VECTOR 
C 

CALL  GAUSSJ  (ALSS,NB,140,RSS.IER) 

C 

C  IF  ERROR  -  BOMB  THE  TURKEY 
C 

IF  (IER.EQ.-l)  GO  TO  9999 
C 

RETURN 

C 

9999  WRITE  (*  ,*)  ’  MATRIX  IS  SINGULAR  -  CHECK  YOUR  INPUT’ 
RETURN 
END 
C 

SUBROUTINE  KERNEL  (Z02.IZ0.IZ) 

C 

Q  ***************** **************** ************************** 

C  THIS  SUBROUTINE  CONTAINS  THE  NON-CAUCHY  SINGULAR 
C  PORTION  OF  THE  GREEN’S  FUNCTION  THAT  IS  USED  TO 
C  SOLVE  THE  PROBLEM 

Q  *********************************************************** 

c 

REALM  M.KAP 

COMPLEX  RKHAT(70),SIHAT(70)21 22.FAC  1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR (70,70) 

COMMON  /GEOM/  R,M,KAP,ALPHA,Z1,Z2,STRESS,PHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70),N2,H,H2,N2P2,FAC1FAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJFLGJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 
COMPLEX  ZETA2ETA02ETAB2ETAOBEYE,A,CEMOZ020BINV2ZO 
COMPLEX  ZOBOM,YEPO,XEPO,T1  ,T1  A, T1  B,T1C,T1  D.T1TIME 
COMPLEX  T2,T3,T4,T5,T6,T6A,T6B,T6C,T6D,T6TIME 
COMPLEX  T7.T7 A,T7B,T7C,T7D,T7TIMEFHI  1  .PH12.PHI3FHI4 
COMPLEX  TRY 
C 

C  CALCULATE  CONSTANTS  TO  BE  USED  IN  THE  CALCULATIONS 
C  STARTING  BY  TRANSFORMING  Z  AND  Z0  INTO  ZETA  AND  ZETAO 
C 

TRY  =  CSQRT(Z**2-4.*R*R*M) 

C 

C  MAKE  SURE  YOU  PICK  THE  RIGHT  ROOT 
C 

IF  (REAL(Z).LE.O..AND.REAL(TRY).GT.O.)  TRY  =  -l.TRY 
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ZETA  =  (Z+TRY)/(2.*R) 

TRY  =  CSQRT(Z0**2-cmplx(4.*R*R*M,0.)) 

IF  (REAL(ZO)JLE.O..AND.REAL(TRY).GT.O.)  TRY  =  -1  *TRY 
ZETAO  =  (Z0+TRY)/(2.*R) 

NOW  THE  REST  OF  THE  CONSTANTS 

ZETAB  =  CONJG(ZETA) 

ZETAOB  =  CONJG(ZET  AO) 

ZOBINV  =  17ZETA0B 
EMOZO  =  M/ZETAO 
ZOBOM  =  ZETAOB/M 
EYE  =  (0..1.) 

A  =  (ZOB INV -EMOZO)*  (ZOB INV  -ZET  A0)/(Z0B  INV -ZOBOM) 

C  =  (EMOZO-ZOBINV)*(EMOZO-ZOBOM)/(EMOZO-ZETAO) 

YEPO  =  EMOZO-KAP/ZETAOB+M*M/ZETAO-M/(KAP*ZETAOB)-M*C+A/KAP 
XEPO  =EYE*EMOZO-EYE*KAP/ZETAOB-EYE*M*M/ZETAO+ 

!  EYE*M/(KAP*ZETAOB)+EYE*M*C+EYE*A/KAP 
RYEPO  =  REAL(EYE*YEPO) 

RXEPO  =  REAL  (EYE*  XEPO) 

C 

C  CALCULATE  THE  TERMS  IN  THE  KERNEL 
C  I  REALIZE  THAT  THIS  IS  A  MESS,  BUT  YOU  SHOULD  SEE  THE 

C  PAGES  OF  ALGEBRA  THAT  MAKE  THIS  THING  UP.  MIKE  SANT  ARE 

C  CAN  HELP  YOU  WITH  THIS  IF  YOU  NEED  IT 
C 

T1A  =  (-2.*ZET  A*  *  3+ZETA*  *2/ZETA0B+M/ZET  AOB) 

!  /(KAP*  R*ZET  AOB*  (ZET  A-ZOBIN  V)*  *2) 

TIB  =  -(-2  *M*ZETA**3+M*M*ZETA**2/ZETAO+M**3/ZETAO) 

!  /((ZET  A-EMOZO)**2*R*ZET  AO) 

TIC  =  -2  *A*(-ZETA**4+M*ZETA/ZETA0B) 

!  /(R*M*KAP*(ZETA-Z0BINV)**3) 

T1D  =  2  *EYE*M*ZETA/(KAP*R*(1.+M*M/KAP)) 

T1TIME  =  (ZETA**2*(ZETAB**2+M))/(ZETAB*(ZETA**2-M)**3) 

C 

T2  =  KAP*ZET  A/(R*ZETA0B  *  (ZET  A-ZOB  INV)*  (ZET  A*  *2-M)) 

T3  =  -M*ZETA/(R*ZETAO*(ZETA-EMOZO)*(ZETA**2-M)) 

T4  =  M*C*ZETA**2/(R*(ZETA-EMOZO)**2*(ZETA**2-M)) 

T5  =  -EYE/(R*(ZETA**2-M)*(1.+M*M/KAP)) 

C 

T6A  =  1./(KAP*R*ZETA0B*(ZETA-Z0BINV)*ZETA) 

T6B  =  -M/(R*ZETAO*(ZETA-EMOZO)*ZETA) 

T6C  =  -A/(R*M*KAP*(ZETA-Z0BINV)**2) 

T6D  =  -EYE*M/(R*KAP*(1.+M*M/KAP)*ZETA**2) 

T6TTME  =  -(M*ZETA**6-ZETA**4*(3.*M*M+1.)-M*ZETA**2) 

!  /((ZETA**2-M)**3) 

C 

T7A  =  -(2  *ZETA-ZOBINV)/(KAP*ZETAOB*ZETA**2*(ZETA-ZOBINV)**2) 
T7B  =  M*(2.*ZETA-EMOZO)/(ZETAO*ZETA**2*(ZETA-EMOZO)**2) 

T7C  =  2  *A/(M*KAP*(ZETA-Z0BINV)**3) 
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T7D  =  2.*EYE*M/(KAP*ZETA**3*(1.+M**2/KAP)) 

T7TIME  =  -ZET  A’  *  3*  ( 1  .+M*ZET  A*  *  2)/(R*  (ZET  A*  *2-M)**  2) 

C 

PHI1  =  T2/KAP**2 
PHI2  =  T3 

Pffl3  =  -A*ZETA**2/(R*M*KAP*(ZETA-Z0BINV)**2*(ZETA**2-M)) 

PHI4  =  T5*M/KAP 
C 

C  NOW  TO  CALCULATE  THE  NORMAL  AND  SHEAR  COMPONENTS  OF 
C  THE  KERNELS.  THE  NORMAL  COMPONENT  IS  THE  REAL  PART 

C  OF  PHLPRIME.  THE  "SHEAR"  COMPONENT  IS  THE  OTHER 

C  PART  OF  MUSHKELISHVILI’S  POTENTIAL  EQUATION 
C 

T1  =  (T1A+T1B+T1C+T1D*RYEP0)*T1TIME 
T6  =  (T6A+T6B+T6C+T6D*RYEP0)*T6TIME 
T7  =  (T7A+T7B+T7C+T7D*RYEP0)*T7TIME 
YSHEAR(IZO,IZ)  =  T 1  +T2+T3+T4+T5*R  YEPOfT6+T7 
C 

T1  =  (-EYE*T1A-EYE*T1B+EYE*T1C+T1D*RXEP0)*T1TIME 
T6  =  (-EYE*T6A-EYE*T6B+EYE*T6C+T6D*RXEP0)*T6TIME 
T7  =  (-EYE*T7A-EYE*T7B+EYE*T7C+T7D*RXEPO)*T7TIME 
XSHEAR(IZO,IZ)  =  Tl+EYE*T2+EYE*T3-EYE*T4+T5*RXEPO+T6+T7 
C 

YNORM(IZO.IZ)  =  2.*REAL(PHI1+PHI2+PHI3+PHI4*RYEP0) 

XNORM(IZO,IZ)  =  2.*REAL(-EYE*PHI1-EYE*PHI2+EYE*PHI3+PHI4*RXEP0) 
C 

RETURN 

END 


COMPLEX  FUNCTION  PFORCE(Z) 

C 

Q  ************************************************************ 

C  SUBROUTINE  THAT  CALCULATES  THE  PHI  POTENTIAL  FOR  THE  FORCE 
C  DUE  TO  THE  ELLIPSE,  TO  BE  USED  IN  BUILD  TO  GENERATE  THE 
C  RIGHT  HAND  SIDE  OF  THE  EWUATION.  SEE  PAGE  349  IN 
C  MUSHKHELISHVILI’S  BOOK 

£  ************************************************************ 

c 

REAL*4  M,KAP 

COMPLEX  RKHAT(70),SIHAT(70)Z:i Z2.FAC  1.FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM 1  R,M,KAP,ALPHA^1^2,STRESSEHIIN 
COMMON  /GEOM2/  ITIP,THETA,DZ,PHI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(I40),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70)N2,H,H2,N2P2,FAC1EAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6/  EJ,FJ,GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 
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COMPLEX  ZETA£YE,Z,DPHI 
COMPLEX  TRY 


GET  ZETA 

TRY  =  CSQRT(Z**2-cmplx(4.*R*R*M,0.)) 

IF  (REAL(Z).LE.O..AND.REAL(TRY).GT.O.)  TRY  =  -l.*TRY 
ZETA  =  (Z+TRY)/(2.*R) 

EYE  =  CMPLX(0.,1.) 

CALCULATE  THE  FORCE  COMPONENTS  FROM  THE  POTENTIALS 

EPS  =  M*SIN(2.*ALPHA)*(1.+KAP)/(4.*(M*M+KAP)) 

DPHI  =  ZETA**2/(4.*(ZETA**2-M))-(2  *M*EPS*EYE+M/4.- 
!  .5*(COS(2*  ALPHA)+EYE*SIN(2  *  ALPHA)))/(KAP*(ZETA*  *2-M)) 

PFORCE  =  cmplx(2.*reaI(DPHI),0.) 

RETURN 

END 

COMPLEX  FUNCTION  SFORCE(Z) 

************************************************************ 

CALCULATES  THE  SI  POTENTIAL  FROM  PAGE  349  IN  MUSHKHELISHVILI 

**********************************************************«* 


REALM  M.KAP 

COMPLEX  RKHAT(70),SIHAT(70)Z1  JZ2.FAC  1  ,FAC2 
COMPLEX  XSHEAR(70,70),YSHEAR(70,70) 

COMMON  /GEOM /  R>LKAP,ALPHA,Z1,Z2,STRESS,PHIIN 
COMMON  /GEOM2/  ITIP, THETA, DZJ?HI(70) 

COMMON  /MATBLK/  ALHS(140,140),RHS(140),RSS(140),ALSS(140,140) 
COMMON  /BLOK3/  SI(70),RK(70),N2,H,H2,N2P2,FAC1.FAC2,ABSFAC 
COMMON  /BLOK4/  PSI(70),TAU(70) 

COMMON  /BLOK5/  VI(70),RKHAT,WIK(70,70),SIHAT 
COMMON  /BLOK6 /  EJ.FJ.GJ 

COMMON  /KERNEL/  XNORM(70,70),YNORM(70,70),XSHEAR,YSHEAR 
COMPLEX  ZETA£YE^BPHI,SIPRIM,PHITIM,T1  ,T2,T3,T4 
COMPLEX  T5,T6,T7Z:ZETAB2SMMA 
COMPLEX  TRY 

GET  ZETA 

TRY  =  CSQRT(Z**2-4.*R*R*M) 

IF  (REAL(Z).LE.O..AND.REAL(TRY),GT.O.)  TRY  =  -l.*TRY 
ZETA  =  (Z+TRY)/(2.*R) 

EYE  =  CMPLX(0.,1.) 

ZETAB  =  CONJG(ZETA) 

CALCULATE  THE  FORCE  COMPONENTS  FROM  THE  POTENTIALS 
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EPS  =  M*SIN(2.*ALPHA)*(1.+KAP)/(4.*(M*M+KAP)) 

SIN  ALP  =  SIN(2.*ALPHA) 

COSALP  =  COS(2.*ALPHA) 

ZSMM  =  ZETA**2-M 

PHTTIM  =  ZETA**3*(ZETAB**2+M)/(ZETAB*ZSMM**2) 

A  =  2  *M*EPS*EYE+M/4.-(COSALP+EYE*SINALP)/2. 

ZBPHI  =  PHITIM*(.5-ZETA**2/(2.*ZSMM)+2.*A/(KAP*ZSMM)) 
T1  =  -ZETA**2*(COSALP-EYE*SINALP)/(2.*ZSMM) 

T2  =  -2  *EPS*EYE/ZSMM 
T3  =  -KAP/(4.*ZSMM) 

T4  =  -ZETA**2*(l.+M**2)/(4  *ZSMM**2) 

T5  =  ZETA**4*(1.+M**2)/(2.*ZSMM**3) 

T6  =  2.*A*M*ZETA**2/(KAP*ZSMM**2) 

T7  =  -A*(1.+M*ZETA**2)*(3.*ZETA**2-M)/(KAP*ZSMM**3) 
SIPRIM  =  T1+T2+T3+T4+T5+T6+T7 
SFORCE  =  ZBPHI+SIPRIM 

RETURN 

END 


SUBROUTINE  GAUSSJ  (A,N,NP,B,IER) 

C 

Q  t ********************************************************* 

C  GAUSS-JORDAN  MATRIX  INVERSION  WITH  FULL  PIVOTING  EXCERPTED 
C  FROM  "NUMERICAL  RECIPES  IN  FORTRAN" 

£  ********************************************************** 

c 

DIMENSION  A(NP,NP),B(NP),IPIV(140),INDXR(140),INDXC(140) 

DO  11  J=1,N 
IPIV(J)  =  0 

11  CONTINUE 
DO  22  1=1, N 

BIG  =  0. 

DO  13  J  =  INI 

IF  (IPIV(J).NE.l)  THEN 
DO  12  K=1,N 

IF  (IPIV(K).EQ.O)  THEN 

IF  (ABS(A(J,K)).GE.BIG)  THEN 
BIG  =  ABS(A(J,K)) 

IROW  =  J 
ICOL  =  K 
ENDIF 

ELSEIF  (IPIV(K).GT.  1)  THEN 

WRITE  (*,*)  *  SINGULAR  MATRIX’ 

IER  =  -1 
GO  TO  999 
ENDIF 

12  CONTINUE 
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END  IF 

13  CONTINUE 
IPIV(ICOL)  =  IPIV(ICOL)+l 
IF  (IROW.NE.ICOL)  THEN 

DO  14  L=1N 

DUM  =  A(IROWU) 

A(IROWU)  =  A(ICOLU) 

A(ICOLU)  =  DUM 

14  CONTINUE 
DUM  =  BflROW) 

B(IROW)  =  B(ICOL) 

B(ICOL)  =  DUM 

ENDIF 

INDXR(I)  =  IROW 
INDXC(I)  =  ICOL 
IF  (A(ICOLJCOL)JEQ.O.)  THEN 

WRITE  (*  *)  ’SINGULAR  MATRIX’ 

IER  =  -1 
GO  TO  999 
ENDIF 

PIVINV  =  17A(ICOL,ICOL) 

A(ICOL.ICOL)  =  1. 

DO  16  L  =  1J9 

A(ICOL,L)  =  A(ICOLL)*PIVINV 
16  CONTINUE 

B(ICOL)  =  B(ICOL)*PIVINV 
DO  21  LL=1,N 

IF  (LL.NE.ICOL)  THEN 
DUM  =  A(LL,ICOL) 

A(LL,ICOL)  =  0. 

DO  18  L=1  J9 

A(LL,L)  =  A(LL,L)-A(ICOLL)*DUM 
18  CONTINUE 

B(LL)  =  B(LL)-B(ICOL)*DUM 
ENDIF 

21  CONTINUE 

22  CONTINUE 
DO  24  L  =  N.1,-1 

IF  (INDXR(L).NE.INDXC(L))  THEN 
DO  23  K=1,N 

DUM  =  A(K,INDXR(L)) 

A(K,INDXR(L))  =  A(K,INDXC(L)) 
A(K,INDXC(L))  =  DUM 

23  CONTINUE 
ENDIF 

24  CONTINUE 
999  CONTINUE 

RETURN 

END 
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SUBROUTINE  RUSSIAN (N2P2) 

C 

Q  *********************************************************** 

C  SUBROUTINE  THAT  CALCULATES  THE  POTENTIALS  FOR  A  CURVED  CRACK 
C  FROM  A  FIRST  ORDER  PERTURBATION  SOLUTION.  GIVEN  IN  A  PAPER 
C  BY  GOLDSTEIN  AND  SALAGANIK,  IN  RUSSIAN,  IN  THE  IZVESTIA  MTT, 

C  DM  1970.  SEE  MY  THESIS  TO  GET  THE  FULL  CITATION 

Q  ***************************** ******* *********************** 

c 

COMMON  /SPUN/  X(70), Y(70), Y 1  (70), Y2(70) 

COMMON  /RUSS/  YSTAR(70,70),YPSTAR(70,70),CAPY(70,70) 

C 

DO  10  IT  =  1.N2P2-1 
DO  10  ITAU  =  1.N2P2-1 
IF  (IT.EQ.IT AU)  THEN 
YSTAR(IT,ITAU)  =  Yl(IT) 

YPSTAR(IT,ITAU)  =  Y2(IT) 

C  APY  (IT.ITAU)  =  0. 

ELSE 

Y  ST  AR(IT,IT  AU)  =  (Y(IT)-Y(ITAU))/(X(IT)-X(ITAU)) 

YPSTAR(IT,ITAU)  =  (Y1(IT)-Y1(ITAU))/(X(IT)-X(ITAU)) 

CAPY(IT,ITAU)  =  (Y 1  (IT)- Y  STAR(IT,IT  AU)) 

!  /(X(IT)-X(ITAU)) 

ENDIF 

10  CONTINUE 
RETURN 
END 
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1  Director 

Benet  Weapons  Laboratory  1 

U.S.  Army  Armament  Research, 

Development,  and  Engineering  Center 
ATTN:  SMCAR-CCB-TL 
Watervliet,  NY  12189-4050 


Commander 

U.S.  Army  Missile  Command 
ATTN:  AMSMI-RD-CS-R  (DOC) 

Redstone  Arsenal,  AL  35898-5010 

Commander 

U.S.  Army  Tank- Automotive  Command 
ATTN:  ASQNC-TAC-DIT  (Technical 
Information  Center) 

Warren,  Ml  48397-5000 

Director 

U.S.  Army  TRADOC  Analysis  Command 
ATTN:  ATRC-WSR 

White  Sands  Missile  Range,  NM  88002-5502 
Commandant 

U.S.  Army  Field  Artillery  School 

ATTN:  ATSF-CSI 

Ft.  Sill,  OK  73503-5000 

Commandant 

U.S.  Army  Infantry  School 

ATTN:  ATSH-CD  (Security  Mgr.) 

Fort  Benning,  GA  31905-5660 

Commandant 
U.S.  Army  Infantry  School 
ATTN:  ATSH-CD-CSO-OR 
Fort  Benning,  GA  31905-5660 

Air  Force  Armament  Laboratory 

ATTN:  WL/MNOI 

Eglin  AFB,  FL  32542-5000 

Aberdeen  Proving  Ground 


(Unciass.  oniy)i  Commander 

U.S.  Army  Armament,  Munitions 
and  Chemical  Command 
ATTN.  AMSMC-IMF-L 
Rock  Island,  IL  61299-5000 


1 


Director 

U.S.  Army  Aviation  Research 
and  Technology  Activity 
ATTN:  SAVRT-R  (Library) 
M/S  219-3 

Ames  Research  Center 
Moffett  Field,  CA  94035-1000 


2  Dir,  USAMSAA 
ATTN:  AMXSY-D 

AMXSY-MP,  H.  Cohen 

1  Cdr,  USATECOM 
ATTN:  AMSTE-TC 

3  Cdr,  CRDEC,  AMCCOM 
ATTN:  SMCCR-RSP-A 

SMCCR-MU 

SMCCR-MSI 

1  Dir,  VLAMO 

ATTN:  AMSLC-VL-D 


10  Dir,  BRL 

ATTN:  SLCBR-DD-T 


119 


No.  Of 

Copies  Organization 

2  PEO-Armaments 

ATTN:  SFAE-AR-PM, 

D.  Adams 

T.  McWilliams 

Picatinny  Arsenal,  NJ  07806-5000 

1  U.S.  Army  Belvoir  RD&E  Center 
ATTN:  STRBE-JBC,  C.  Kominos 
Fort  Belvoir,  VA  22060-5606 

1  Commander 
DARPA 

ATTN:  J.  Kelly 
1400  Wilson  Blvd. 

Arlington,  VA  22209 

2  Materials  Technology  Laboratory 
ATTN:  SLCMT-MEC, 

B.  Halpin 
T.  Chou 

Watertown,  MA  02172-0001 

1  Commander 

U.S.  Army  Laboratory  Command 
ATTN:  AMSLC-TD,  R.  Vitali 
Adelphi,  MD  20783-1145 

1  Commander 

U.S.  Army  Missile  Command 
ATTN.  AMSMI-RD,  W.  McCorkle 
Redstone  Arsenal,  AL  35898-5010 

2  Commander 

U.S.  Army  Laboratory  Command 
Harry  Diamond  Laboratory 
ATTN:  SLCHD-TS-NT,  A.  Frydman 
2800  Powder  Mill  Road 
Adelphi,  MD  20783-1197 


No.  of 

Copies  Organization 
9  Director 

Benet  Weapons  Laboratory 
U.S.  Army  Armament  Research, 
Development,  and  Engineering  Center 
ATTN:  SMCAR-CCB, 

V.  Montuori 
J.  Keane 
T.  Allen 
J.  Vasilakis 
G.  Friar 
J.  Zweig 

L.  Johnson 
T.  Simkins 

J.  Wrzockalski 
Watervliet,  NY  12189-5000 

7  Commander 

U.S.  Army  Armament  Research, 
Development,  and  Engineering  Center 
ATTN:  SMCAR-CCH-T, 

S.  Musalii 
P.  Christian 

K.  Fehsal 

SMCAR-CCH-V,  E.  Fennell 
SMCAR-CCH,  J.  DeLorenzo 
SMCAR-CC, 

R.  Price 
J.  Hedderich 

Picatinny  Arsenal,  NJ  07806-5000 

2  Commander 

U.S.  Army  Armament  Research, 
Development,  and  Engineering  Center 
ATTN:  SMCAR-TD, 

M.  Lindner 

T.  Davidson 

Picatinny  Arsenal,  NJ  07806-5000 

1  Commander 

Production  Base  Mod.  Agency 

U.S.  Army  Armament  Research, 
Development,  and  Engineering  Center 
ATTN:  AMSMC-PBM-K 
Picatinny  Arsenal,  NJ  07806-5000 
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No.  oi 
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3  PEO-Armaments 
Project  Manager 

Tank  Main  Armament  Systems 
ATTN:  SFAE-AR-TMA,  COL  Hartline 
SFAE-AR-TMA-MD,  C.  Kimker 
SFAE-AR-TMA-ME,  K.  Russell 
Picatinny  Arsenal,  NJ  07806-5000 

2  Commander 

Wright-Patterson  Air  Force  Base 
ATTN:  AFWAML, 

J.  Whitney 
R.  Kim 

Dayton,  OH  45433 

2  Pacific  Northwest  Laboratory 

A  Division  of  Battelle  Memorial  Institute 
ATTN:  M.  Smith 
M.  Garnich 
P.O.  Box  999 
Richland,  WA  99352 

2  David  Taylor  Research  Center 
ATTN:  R.  Rockwell 
W.  Phyillaier 

Bethesda,  MD  20054-5000 
1  Director 

Los  Alamos  National  Laboratory 

ATTN:  D.  Rabern 

WX-4  Division,  Mail  Stop  G-787 

P.O.  Box  1663 

Los  Alamos,  NM  87545 

4  Director 

Sandia  National  Laboratories 
Applied  Mechanics  Department, 
Division-8241 
ATTN:  C.  Robinson 
G.  Benedetti 

K.  Perano 
W.  Kawahara 
P.O.  Box  969 

Livermore,  CA  94550-0096 


3  Director 

Lawrence  Livermore  National 
Laboratory 

ATTN:  R.  Christensen 
W. Feng 
S.  deTeresa 
P.O.  Box  808 
Livermore,  CA  94550 

3  University  of  Delaware 

Center  for  Composite  Materials 
ATTN:  J.  Gillespe 
B.  Pipes 
M.  Santare 

201  Spencer  Laboratory 
Newark,  DE  19716 

2  Civil  Engineering  Department 
N.C.  State  University 
P.O.  Box  7908 
ATTN:  W.  Rasdorf 

L.  Spainhour 
Raleigh,  NC  27696-7908 

1  Pennsylvania  State  University 

Department  of  Engineering  Science 
and  Mechanics 
ATTN:  T.  Hahn 
227  Hammond  Building 
University  Park,  PA  16802 

1  University  of  Utah 

Department  of  Mechical  and 
Industrial  Engineering 
ATTN:  S.  Swanson 
Salt  Lake  City,  UT  84112 

1  Stanford  University 

Department  of  Aeronautics  and 
Aeroballistics  Durant  Bldg. 

ATTN:  S.  Tsai 
Stanford,  CA  94305 

1  Custon  Analytical  Engineering 

Systems,  Inc. 

ATTN:  A.  Alexander 
Star  Route  Box  4A 
Flintstone,  MD  21530 


121 


Zak  Technologies,  Inc. 

ATTN:  A.  Zak 
2310  Belmore  Drive 
Champaign,  IL  61820 

Chamberlain  Manufacturing  Corp. 
Waterloo  Facility 
ATTN:  T.  Lynch 
550  Ester  St. 

P.O.  Box  2335 
Waterloo,  IA  50704 

Olin  Corporation 
Flinchbaugh  Division 
ATTN:  E.  Steiner 
B.  Stewart 
P.O.  Box  127 
Red  Lion,  PA  17356 

Olin  Corporation 
ATTN:  L.  Whitmore 
10101  9th  St.  North 
St.  Petersburg,  FL  33702 

Alliant  Techsystems,  Inc. 

ATTN:  C.  Candland 

J.  Bode 

K.  Ward 

5640  Smetana  Drive 
Minnetonka,  MN  55343 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it 
publishes.  Your  comments/answers  below  will  aid  us  in  our  efforts. 

1.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of 
interest  for  which  the  report  will  be  used.) _ 


2.  How,  specifically,  is  the  report  being  used?  (Information  source,  design  data,  procedure, 
source  of  ideas,  etc.) _ 


3.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or 
dollars  saved,  operating  costs  avoided,  or  efficiencies  achieved,  etc?  If  so,  please 
elaborate.  _  _ 


4.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports? 
(Indicate  changes  to  organization,  technical  content,  format,  etc.)  _ 


BRL  Report  Number  BRL-tr-3277 _  Division  Symbol 

Check  here  if  desire  to  be  removed  from  distribution  list.  _ 

Check  here  for  address  change.  _ 

Current  address:  Organization  _ 

Address _ 
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ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  21005-5066 


